Compositions and methods for multiplexed quantitative analysis of cell lineages

ABSTRACT

Compositions and methods are provided for measuring population size for a plurality of clonal cell populations in the same individual, e.g., for measuring tumor size for a plurality of clonally independent tumors within the same individual. A subject method can include: (a) contacting an individual with a plurality of cell markers that are heritable and distinguishable from one another, to generate a plurality of distinguishable lineages of heritably marked cells; (b) after sufficient time has passed for the heritably marked cells to undergo at least one round of division, detecting and measuring quantities of at least two of the plurality of cell markers present in the contacted tissue, thereby generating a set of measured values; and (c) using the set of measured values to calculate the number of heritably marked cells that are present (e.g., for at least two of the distinguishable lineages of heritably marked cells).

CROSS-REFERENCE

This application is a continuation of U.S. patent application Ser. No.15/940,818 filed Mar. 29, 2018, which claims the benefit of U.S.Provisional Patent Application No. 62/481,067 filed Apr. 3, 2017, whichapplications are incorporated herein by reference in their entirety.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

This invention was made with Government support under contractsCA124435, CA194910, CA207133 and GM118165 awarded by the NationalInstitutes of Health. The Government has certain rights in theinvention.

INTRODUCTION

Genome sequencing has catalogued the somatic alterations in humancancers at the genome-wide level and identified many potentiallyimportant genes (e.g., putative tumor suppressor genes, putativeoncogenes, genes that could lead to treatment resistance orsensitivity). However, the identification of genomic alterations doesnot necessarily indicate their functional importance in cancer, and theimpact of gene inactivation or alteration, alone or in combination withother genetic alterations (either somatic or germline) ormicroenvironmental differences, remains difficult to glean from cancergenome sequencing data alone.

The molecular and cellular impacts of genetic alterations on neoplasticgrowth have been directly investigated using knockdown, knockout, andoverexpression studies in cell lines as well as genetically engineeredmouse model systems. Over the past several decades the analysis of genefunction in cancer cell lines in culture has provided insights into manyaspects of cancer. However, the near-optimal growth of cancer cell linesin culture, widespread pre-existing genetic and epigenetic changes, andthe lack of the autochthonous microenvironment limit the ability ofthese systems to provide insight into how different genes constrain ordrive in vivo phenotypes (e.g. cancer growth, metastasis, therapyresponses). In contrast, genetically engineered mouse models of humancancer facilitate the introduction of defined genetic alterations intonormal adult cells which results in the initiation and growth of tumorswithin their natural in vivo setting. This is of particular importanceas many pathways are influenced by properties of the in vivo tumormicroenvironment.

While in vivo systems such as CRISPR/Cas based genetic targeting haveincreased the scale of in vitro and in vivo functional analyses, in vivosystems have continued to rely on relatively crude measurements of tumorgrowth, limiting their application to the analysis of genes with themost dramatic effects. The lack of rigorously quantitative systems toanalyze gene function in vivo has precluded a broad understanding ofpathways that drive or constrain tumor growth, or impact any of theother important aspects of carcinogenesis (e.g., tumor suppressorpathways).

There is a need for compositions and methods that facilitate precisequantification of clonal population size (e.g., the size of each tumor,the number of neoplastic cells in each tumor or subclone, and the like)in an individual with a plurality of clonal cell populations (e.g., aplurality of distinguishable cell lineages—being either distinct,identifiable tumors, or distinct identifiable subclones within a tumor).The compositions and methods of this disclosure address this need, andprovide the ability to uncover whether different individual genes (e.g.,tumor suppressors, oncogenes) or genetic alterations (e.g. insertions,deletions, point mutations), or combinations of genes and/or geneticalterations, have different overall effects on cell population growth(e.g., tumor growth), as well as other phenotypes of importance (e.g.,tumor evolution, progression, metastatic proclivity). The compositionsand methods of this disclosure also provide the ability to test theeffect of potential therapeutics, e.g., radiation, chemotherapy,fasting, compounds such as drugs, biologics, etc., on the growth ofmultiple different clonal cell populations (e.g., multiple tumors ofsimilar genotype but with different initiation events, multiple tumorsthat have different genotypes, and the like) within the same tissue(e.g., within the same individual), which would drastically reduce errorintroduced by sample-to-sample variability (e.g., animal-to-animalvariability). These methods also facilitate development and testing ofrational drug combinations.

SUMMARY

Compositions and methods are provided for measuring population size fora plurality of clonal cell populations in the same tissue (e.g., in thesame individual) or in different tissues. As an example, in some cases asubject method is a method of measuring tumor size for a plurality ofclonally independent tumor cell populations (e.g., different tumors) inthe same tissue (e.g., in the same individual).

As an illustrative example, as described below in the working examples,the inventors combined cell barcoding (e.g., tumor barcoding) andhigh-throughput sequencing (referred to in the working examples as“Tuba-seq”) with genetically engineered mouse models of human cancer toquantify tumor growth with unprecedented resolution. Precisequantification of individual tumor sizes allowed them to uncover theimpact of inactivating different tumor suppressor genes (e.g., knowntumor suppressor genes). Further, the inventors integrated these methodswith multiplexed CRISPR/Cas9-mediated genome editing, which allowedparallel inactivation and functional quantification of a panel ofputative tumor suppressor genes—and led to the identification offunctional lung tumor suppressors. The method is a rapid, multiplexed,and highly quantitative platform to study the impact of geneticalterations on cancer growth in vivo.

Also as described in the working examples below, the inventors usedmultiplexed somatic homology directed repair (HDR) with barcoded HDRdonor templates to produce genetically diverse barcoded tumors (e.g.,tumors that have genetically diverse point mutations in a defined gene)within individual mice, and employed quantitative tumor analysis (usinghigh-throughput sequencing) to rapidly and quantitatively interrogatethe function of multiple precise mutations (e.g., defined pointmutations) simultaneously in the same animal.

In some embodiments, a subject method includes a step of contacting atissue (e.g., muscle, lung, bronchus, pancreas, breast, liver, bileduct, gallbladder, kidney, spleen, blood, gut, brain, bone, bladder,prostate, ovary, eye, nose, tongue, mouth, pharynx, larynx, thyroid,fat, esophagus, stomach, small intestine, colon, rectum, adrenal gland,soft tissue, smooth muscle, vasculature, cartilage, lymphatics,prostate, heart, skin, retina, and the reproductive and genital systems,e.g., testicle, reproductive tissue, etc.) with a plurality of cellmarkers that are heritable and distinguishable from one another, togenerate a plurality of distinguishable lineages of heritably markedcells within the contacted tissue. In some embodiments, the cell markersused to contact the tissue are barcoded nucleic acids (e.g., RNAmolecules; or circular or linear DNA molecules such as plasmids, naturalor synthesized single- or double-stranded nucleic acid fragments, andminicircles). In some embodiments (e.g., in cases where the cell markersare barcoded nucleic acids), the cell markers can be delivered to thetissue via viral vectors (e.g., lentiviral vectors, adenoviral vectors,adeno-associated viral (AAV) vectors, and retroviral vectors). In somecases, the tissue to be contacted already includes neoplastic cellsprior to contact with cell markers. In some cases, the cell markers caninduce neoplastic cell formation and/or tumor formation. In some cases,components linked to the cell markers can induce neoplastic cellformation and/or tumor formation. In some cases, the cell markers arebarcoded nucleic acids that can induce neoplastic cell formation and/ortumor formation (e.g., homology directed repair (HDR) DNA donortemplates; nucleic acids encoding a genome editing protein(s); nucleicacids encoding oncogenes; nucleic acids encoding a protein(s), e.g.,wild type and/or mutant protein(s) [e.g., wild type or mutant cDNA thatencodes a protein that is detrimental to tumors, e.g., in some way otherthan growth/proliferation]; CRISPR/Cas guide RNAs; short hairpin RNAs(shRNAs); nucleic acids encoding targeting components for other genomeediting systems; etc.).

Subject methods can also include (after sufficient time has passed forat least a portion of the heritably marked cells to undergo at least oneround of division) a step of detecting and measuring quantities of atleast two of the plurality of cell markers present in the contactedtissue—thereby generating a set of measured values, which represent theidentity and quantity of cell markers that remain in the contactedtissue, e.g., heritably associated with the marked cells. In some cases(e.g., when the cell markers are barcoded nucleic acids) the detectingand measuring can be performed via a method that includeshigh-throughput sequencing and quantification of the number of sequencereads for each detected barcode.

In some cases, the generated set of measured values is used as input tocalculate (e.g., using a computer) the number of heritably marked cellspresent in the contacted tissue (e.g., for at least 2, at least 3, atleast 4, at least 5, at least 100, at least 1,000, at least 10,000, orat least 100,000 of the detected distinguishable lineages of heritablymarked cells)(e.g., in some cases in a range of from 10 to 1,000,000;from 10 to 100,000; from 10 to 10,000; or from 10 to 1,000; of thedetected distinguishable lineages of heritably marked cells). Thecalculated number of heritably marked cells can be absolute (e.g., anactual number of cells determined to be present), or can be relative(e.g., a population size for a first lineage of heritably marked cellscan be determined relative to a population size for a second lineage ofheritably marked cells without necessarily determining the actual numberof cells present in either lineage).

In some embodiments, a subject method includes a step of administering atest compound (e.g., a drug) to the tissue (e.g., via administration toan individual, via contacting a synthetic ex vivo tissue such as anorganoid, and the like), e.g., after introducing the cell markers, e.g.,after a step of inducing neoplastic cells (or subclones) via contactingtissue with the plurality of cell markers. In some such cases, the stepof administering the test compound is followed by a step of measuringpopulation size (e.g., tumor size, number of neoplastic cells in eachtumor) for a plurality of marked cell lineages/cell populations. Becausemultiple cell populations can be measured (e.g., multiple tumor sizescan be measured) for distinct and distinguishable marked cell lineageswithin the same tissue (e.g. within the same animal), the risk of errordue to sample-to-sample variation (e.g., animal-to-animal variation) ofdrug response can be greatly reduced, if not eliminated.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention is best understood from the following detailed descriptionwhen read in conjunction with the accompanying drawings. The patent orapplication file contains at least one drawing executed in color. Copiesof this patent or patent application publication with color drawing(s)will be provided by the Office upon request and payment of the necessaryfee. It is emphasized that, according to common practice, the variousfeatures of the drawings are not to-scale. On the contrary, thedimensions of the various features are arbitrarily expanded or reducedfor clarity. Included in the drawings are the following figures.

FIG. 1. Tuba-seq combines tumor barcoding with high-throughputsequencing to allow parallel quantification of tumor sizes. a, Schematicof Tuba-seq pipeline to assess lung tumor size distributions. Tumorswere initiated in KrasLSL-G12D/+; Rosa26LSL-Tomato (KT), KT;Lkb1flox/flox (KLT), and KT; p53flox/flox (KPT) mice with Lenti-mBC/Cre,a virus containing a random 15-nucleotide DNA barcode (BC). Tumor sizeswere calculated via bulk barcode sequencing of the DNA from the tumorbearing lungs. b, Fluorescence dissecting scope images of lung lobesfrom KT, KLT, and KPT mice with Lenti-mBC/Cre initiated tumors. Lunglobes are outlined with white dashed lines. The titer of Lenti-mBC/Creis indicated. Different titers were used in different genetic backgroundto generate approximately equal total tumor burden despite differencesin overall tumor growth. Scale bars in upper panels=5 mm. Scale bars inlower panels=1 mm. c, Tumor size distributions in KT, KLT, and KPT mice(number of mice per group is indicated). Each dot represents a tumor.The area of each dot is proportional to the number of cancer cells ineach tumor. A dot corresponding to the approximate number of cancercells in a spherical tumor with a 1 mm diameter is shown to the right ofthe data for reference.

FIG. 2. Tuba-seq is a robust and reproducible method to quantify tumorsizes. a, DADA2, a denoising algorithm designed for deep sequencing ofamplicon data, eliminates recurrent read errors that can appear asspurious tumors. Cell lines with known barcodes were added to each lungsample from each mouse (5×10⁵ cells each). Recurrent read errors thatderive from these known barcodes appear as spurious tumors at ˜5,000cells. DADA2 identifies and greatly reduces these recurrent read(sequencing) errors. b,c, Technical replicate sequencing librariesprepared from an individual bulk lung sample demonstrate highcorrespondence between individual lesion sizes (b) and size profiles (c)(tumors at the 50 to 99.9th percentiles are shown). d, Our analysispipeline is robust to variation in read depth, GC content of the DNAbarcodes, and diversity of the barcode library. Tumors were partitionedinto thirds corresponding to high, moderate, and low levels of eachtechnical parameter: the sequencing depth, GC content of tumor barcodes,and estimated number of unique barcodes (see Methods). Whiskers cappedat 1.5 IQR. e, Reproducibility of size distributions across five KLTmice. Mice have overall similar size profiles despite smallmouse-to-mouse differences in tumor sizes. Sizes of the tumors at theindicated percentiles in individual mice are connected by a line. f,Reproducibility of size profiles improves when tumors within the samemouse are compared, suggesting significant mouse-to-mouse variability intumor sizes. Tumors in each mouse were partitioned into two groups andthe profiles of these groups were compared. Sizes of the tumors at theindicated percentiles in an individual mouse are connected by a line. g,Unsupervised hierarchical clustering of the KT, KPT, and KLT mice basedon the total least-squares distance between tumors sizes at definedpercentiles (clustered by Ward's Variance Minimization Algorithm). Micecluster by genotype suggesting that Tuba-seq identifies reproducibledifferences in the size spectrum of each genotype.

FIG. 3. Massively parallel quantification of tumor sizes enablesprobability distribution fitting across multiple genotypes. a, b, Tumorsize at the indicated percentile in KLT (n=5) mice (a) and KPT (n=3)mice (b) versus tumor size at the indicated percentile in KT mice (n=7).Each percentile was calculated using all tumors from all mice of eachgenotype 11 weeks after tumor initiation with Lenti-mBC/Cre. c, Tumorsizes at the indicated percentiles for each genotype relative to KTtumors at the same percentiles. Error bars are 95% confidence intervalsobtained via bootstrapping. Percentiles that are significantlydifferently from the corresponding KT percentiles are in color. d, Asanticipated for exponential tumor growth with normally distributedgrowth rates, tumor size distributions were most closely fit by alognormal distribution. Tumors in KLT mice are best described by alognormal distribution throughout their entire size spectrum (middle).The tumor size distributions in KT mice (left) and KPT mice (right) werebetter explained by combining a lognormal distribution at smaller scaleswith a power-law distribution at larger scales. These differences arefundamentally important in considering how individual genes (orcombinations of genes) lead to increased tumors growth. Power-lawrelationships decline linearly on log-log axes, consistent with rare,yet very large tumors within the top ˜1% of tumors in KT mice and ˜10%of tumors in KPT mice. Note: only tumors in KPT mice ever exceed onemillion cells after 11 weeks, consistent with p53-deficiency enablingthe generation of the largest tumors in our study.

FIG. 4. Rapid quantification of tumor suppressor phenotypes usingTuba-seq and multiplexed CRISPR/Cas9 mediated gene inactivation. a,Schematic of the Lenti-sg TS-Pool/Cre vector that contain atwo-component barcode with an 8-nucleotide “sgID” sequence linked toeach sgRNA as well as a random 15 nucleotide random barcode (BC). b,Lenti-sg TS-Pool/Cre contains four vectors with inert sgRNAs and elevenvectors targeting known and candidate tumor suppressor genes. Each sgRNAvector contains a unique sgID and a random barcode. NT=Non-Targeting. c,Schematic of multiplexed CRISPR/Cas9-mediated tumor suppressorinactivation coupled with Tuba-seq to assess the function of eachtargeted gene on lung tumor growth in vivo. Tumors were initiated withLenti-sgTS-Pool/Cre virus in KT and KT; H11^(LSL-Cas9) (KT; Cas9) mice.d, Bright field (top) and fluorescence dissecting scope images (bottom)of lung lobes from KT and KT; Cas9 mice 12 weeks after tumor initiationwith Lenti-sg TS-Pool/Cre. Lung lobes are outlined with white dashedlines in the fluorescence images. Viral titer is indicated. Scale bars=5mm. e, Histology confirms that KT mice have hyperplasias and smalltumors, while KT; Cas9 mice have much larger tumors. Viral titer isindicated. Top scale bars=3 mm. Bottom scales bars=500 μm.

FIG. 5. Tuba-seq uncovers known and novel tumor suppressors withunprecedented resolution. a, Analysis of the relative tumor sizes in KT;Cas9 mice 12 weeks after tumor initiation with Lenti-sgTS-Pool/Creidentified six tumor growth suppressing genes. Relative size of tumorsat the indicated percentiles represents merged data from 8 mice,normalized to the average size of sg Inert tumors. 95% confidenceintervals are shown. Percentiles that are significantly greater than sgInert are in color. b, Estimates of mean tumor size, assuming alognormal tumor size distribution, identified sgRNAs that significantlyincrease growth in KT; Cas9 mice. Bonferroni-corrected, bootstrappedp-values are shown. p-values <0.05 and their corresponding means arebold. c, Relative size of the 95th percentile tumors (left), lognormal(LN) mean (middle), and lognormal (LN) p-value (right) for tumors witheach sgRNA in KT and KT; Cas9 mice 12 weeks after tumor initiation, andKT; Cas9 mice 15 weeks after tumor initiation. d, Fold change in overallsgID representation in KT; Cas9 mice relative to KT mice (AsgIDRepresentation) identified several sgRNAs that increase inrepresentation, consistent with increased growth of tumors withinactivation of the targeted tumor suppressor genes. ΔsgIDRepresentation is the fold change in percent of reads with each sgID inKT; Cas9 mice versus KT mice, normalized such that ΔsgID Representationfor sg Inert=1. Means and 95% confidence intervals are shown. e,f, Theability to detect tumor suppressive effects is drastically improved byincorporating individual tumor sizes from barcode sequencing compared toonly incorporating sgRNA representation. All current approaches rely onsgRNA representation, which is far inferior to Tuba-seq. The relativesize of the 95th percentile tumor and the lognormal statisticalsignificance determined by Tuba-seq identified more genes as tumorsuppressors than the average fold change in ΔsgID representation andtheir associated p-values (e and f). Error bars in (e) are 95%confidence intervals. Dotted lines in (f) indicate the 0.05 significancethreshold. Dot color corresponds to the sgRNA color in FIG. 4 b.

FIG. 6. Independent methods identify Setd2 as a potent suppressor oflung tumor growth. a, The percent of reads containing indels at thetargeted locus was normalized to the average percent of reads containingindels in 3 independent Neomycin loci. This value is plotted versus thesize of the 95^(th) percentile tumor for each sgRNA for three individualmice. We demonstrate a high frequency of indels in Setd2, Lkb1, and Rb1consistent with selection for on-target sgRNA cutting. Each dotrepresents an sgRNA from a single mouse. sgNeo dots are in black and allother dots are colored according to FIG. 4b . b, Fluorescence dissectingscope images and H&E of lung lobes from KT; Cas9 mice infected(transduced) with Lenti-sgSetd2#1/Cre, Lenti-sgSetd2#2/Cre, orLenti-sgNeo2/Cre analyzed 9 weeks after tumor initiation. Lung lobes areoutlined with white dashed lines in the fluorescence dissecting scopeimages. Upper scale bars=5 mm. Lower scale bars=2 mm. c, Quantificationof percent tumor area by histology shows a significant increase in tumorburden in KT; Cas9 mice infected (transduced) with Lenti-sgSetd2#1/Creor Lenti-sgSetd2#2/Cre compared to KT mice infected (transduced) withthe same virus. Each dot represents a mouse and the bars are the mean. *p-value <0.05. NS=not significant. d, Tumor size at the indicatedpercentile from KT; Cas9 mice with Lenti-sgSetd2#1/Cre initiated tumorsversus Lenti-sgNeo2/Cre initiated tumors (N=4 mice/group). Percentileswere calculated using all tumors from all mice in each group

FIG. 7. Frequency of genomic alterations in human lung adenocarcinomaand description of tumor initiation and barcoding. a, The percent oftumors with potentially inactivating alterations (frameshift ornon-synonymous mutations, or genomic loss) in each tumor suppressor geneis shown for all tumors (All) as well as in tumors with oncogenic KRASmutations (KRAS^(mut)). The number and percent of tumors with oncogenicmutations in KRAS in each dataset is indicated. b, Inhalation ofbarcoded lentiviral-Cre vectors initiate lung tumors in geneticallyengineered mouse models. Importantly, the lentiviral vectors stablyintegrate into the genomes of the transduced cells. The relativeexpansion of each uniquely barcoded cell can be determined byhigh-throughput sequencing-based methods. c, Hemotoxilin and Eosin (H&E)staining of lung tissue sections from Kras^(LSL-G12D/+);R26^(LSL-Tomato) (KT) mice infected (transduced) with 1.7×10⁴ Lenti-Crevirus. These mice develop small expansions of neoplastic cells as wellas larger adenomas. Scale bars=50 μm.

FIG. 8. Tuba-seq pipeline to quantify tumor sizes in vivo. a, Illumina®sequencing of the DNA barcode region of the integrated lentiviralvectors enables precise measurement of lesion sizes. First, reads withpoor Phred quality scores or unexpected sequences were discarded. Next,reads were piled-up into groups with unique barcodes. RecurrentIllumina® sequencing errors were delineated from small lesions usingDADA2, a model of Illumina® sequencing errors initially designed toidentify full read-length deep-sequencing amplicons. Small barcodepileups deemed to be recurrent sequencing errors from the amplifiedbarcode region of large tumors were combined with these larger pileupsby this clustering algorithm. Read pileups were translated into absolutecell number using the benchmark controls. Lastly, a minimum cutoff tocall lesions was established using both sequencing information andabsolute cell number to maximize reproducibility of the pipeline. b,c, Aunique read pileup may not correspond to a unique lesion but ratherarise from recurrent sequencing errors of the barcode from a very largetumor (e.g., much larger tumor). DADA2 was used to merge small readpileups with larger lesions of sufficient size and sequence similarity.The algorithm calculates the sequencing error rates from thenon-degenerate regions of our deep sequenced region (i.e. the region ofthe lentiviral vectors that flank the barcode) (b). The likelihood ofevery transition and transversion (A to C shown) was calculated forevery Illumina® Phred score to generate an error model specific for eachrun (c). The advertised Phred error rates (red) are generally lower thanobserved (black; LOESS regression used for regularization). These errormodels (trained to each Illumina® machine) were then used to determineif smaller read pileups should be bundled into larger pileups withstrong sequence similarity (suggesting that the smaller pileup is arecurrent read error) or left as a separate lesion. d-f, We sequencedour first experimental samples (KT, KLT, and KPT from FIG. 1) on threedifferent Illumina® machines to vet and parameterize DADA2. A soundlesion calling protocol was expected to show (d) strong similarity inthe number of called lesions, (e) good correlation between barcodesizes, and (f) similar mean sizes of each sgID pool across the 3 runs.The three runs naturally varied in sequencing depth (40.1×10⁶, 22.2×10⁶,and 34.9×10⁶ reads after pre-processing) and naturally varied in theirexpected error rate per base (0.85%, 0.95%, and 0.25%) offering usefultechnical perturbations to vet concordance of the method. We found thattruncating lesion sizes at 500 cells and truncating the DADA2 clusteringprobability (omega) at 10-10 (red square) offered a profile of lesionsizes at very small scales, while still minimizing variability in ourtest metrics.

FIG. 9. Benchmark controls allow calculation of the number of cancercells in each tumor within each lung sample. a, Schematic of theprotocol using three benchmark control cell lines with known barcodes.5×10⁵ cells of each cell line was added to each lung sample. DNA wasthen extracted from the lung plus all three benchmark controls, and thebarcodes were PCR amplified and deep sequenced. We then calculated thenumber of cancer cells in each tumor within that lung sample by dividingthe % reads associated with the benchmarks by the % reads observed fromeach tumor (unique barcode) and multiplying by 5×10⁵ to obtain cancercell number. b, Example of two lungs with very different tumor burdens.These benchmark cell lines can be used determine the number of cancercells within individual tumors regardless of overall tumor burden. Itshould also be noted, that the surrounding “normal lung” tissue has noimpact on this calculation as this tissue has no lentiviral integrationand thus will contribute no reads. The generation of a titration ofbenchmark controls (e.g., of 5×10⁵, 5×10⁴, 5×10³, 5×10², or 50 cells)facilitates the resolution of Tuba-seq to be extended to smaller clonalexpansions).

FIG. 10. The DADA2-based tumor calling pipeline is robust andreproducible. a, Tumor sizes exhibited a subtle GC-bias. Residual tumorsize variability was minimized by log-transformation of sizes andnormalization of each tumor by the mean size of each sgRNA in everymouse. Barcodes with intermediate GC-content appear to be PCR-amplifiedmost efficiently. A 4th-order polynomial fit to the residual biascorrected lesion sizes most effectively. This correction was calculatedand applied to all subsequent analyses, which adjusted each lesion sizeby an average of 5%, and reduced the standard deviation of lesions sizesof each sgID in each mouse by only 2.9% relative to the mean—suggestingthat, while measurable, variability introduced by GC-bias was minimal.b, The random barcodes exhibited a high-degree of randomness across theintended nucleotides. c, Number of lesions called per mouse usingTuba-seq. Numbers of tumors above two different cell number cutoffs(1000 and 500) are shown as the average number of tumors per mouse ±thestandard deviation. KT mice were exposed to a high titer (6.8×10⁵) (usedin the main text) and a lower titer (1.7×10⁵; KP^(low)). There was nostatistically significant difference in the number of tumors observedper capsid at either cell cutoff suggesting that barcode diversity isstill not limited above half a million tumors and that small tumors arenot caused by tumor crowding. d, Unsupervised hierarchical clustering ofthe KT, KT^(low), KPT, and KLT mice based on the total least-squaresdistance between tumors sizes at defined percentiles (linkage determinedby Ward's Incremental algorithm.) Mice of the same genotype, butdifferent viral titers, cluster together, suggesting that size profiledifferences are determined primarily by tumor genetics (genootype), notdifferences in viral titer. e, f, Lesion sizes are not dramaticallyaffected by differences in read depth. The barcode region from thetumor-bearing lungs of an individual mouse was sequenced at very highdepth and then randomly down-sampled to typical read depth. (e) Thetumor size distributions of the full (x-axis) and downsampled (y-axis)data sets were very similar, indicating our analysis parameters areunbiased by, and fairly robust to, read depth. (f) The percentilescalculations are also reproducible upon downsampling. g, KT, KLT, andKPT mice with Lenti-mBC/Cre initiated tumors (from FIG. 1) have tumorswith six unique Lenti-sg/D-BC/Cre viruses (each harboring a unique sgIDand naturally varying barcode diversity). This allowed us to quantifythe variation in DADA2-called tumor sizes with six replicates withineach mouse. Tumor size distributions are reproducibly called when usingall tumors from each mouse and when using each subset of tumors with agiven sgID. The size of the tumor at the indicated percentiles areplotted for KT (left), KLT (middle), and KPT (right). Each dotrepresents the value of a percentile calculated using tumors within asingle sgID. Percentiles are represented in grey-scale. The sixreplicate percentile values of tumor size with differing sgIDs aredifficult to distinguish since their strong correlation means thatmarkers for each sgID are highly overlapping.

FIG. 11. Efficient genome editing in lung tumors initiated withLentiviral-sgRNA/Cre vectors in mice with the H11^(LSL-Cas9) allele. a,Schematic of the experiment to test somatic genome editing in the lungcancer model using a Lenti-sg Tomato/Cre (Lenti-sgTom/Cre) viral vectorand the H11^(LSL-Cas9) allele. All mice were homozygous for theR26^(LSL-Tomato) allele to determine the frequency of homozygousdeletion. b, Fluorescence dissecting scope images of a lung lobe from aKPT; Cas9 mouse with Lenti-sgTomato/Cre-initiated tumors.Tomato-negative tumors are outlined with dashed lines. Top scale bars=5mm; bottom scale bars=1 mm. c, Immunohistochemistry for Tomato proteinuncovered Tomato-positive (Pos), Tomato-mixed (Mixed), andTomato-negative (Neg) tumors. Tumors are outlined with dashed lines.Scale bars=200 μm. d, Quantification of Tomato expression in four KPT;Cas9 mice with Lenti-sg Tom/Cre initiated tumors indicates thatapproximately half of the tumors have CRISPR/Cas9-mediated homozygousinactivation of the targeted gene in at least a fraction of the cancercells. Percent of Tomato positive, mixed, and negative tumors is shownwith the number of tumors in each group indicated in brackets. e,Schematic of the experiment to test somatic genome editing in the lungusing Lenti-sg Lkb1/Cre virus and the H11^(LSL-Cas9) allele. f,Fluorescence dissecting scope images of lung lobes of KT and KT; Cas9mice infected (transduced) with Lenti-sg Lkb1/Cre show increased tumorburden in the KT; Cas9 mouse. Lung lobes are outlined with white dashedlines. Scale bars=2 mm. g, Tumor burden, represented by lung weight, isincreased in Lenti-sg Lkb1/Cre-infected (transduced) KT; Cas9 micerelative to KT mice, consistent with successful deletion of the tumorsuppressor Lkb1. Normal lung weight is indicated by the dotted redline. * p-value <0.02. Each dot is a mouse and the bar represents themean. h, Western blot showing that Lenti-sg Lkb1/Cre initiated tumors inKT; Cas9 mice express Cas9 and lack Lkb1 protein. Hsp90 shows loading.

FIG. 12. Selection and characterization of sgRNAs targeting eleven knownand candidate tumor suppressor genes. a, sgRNAs were selected based ontheir location within each gene, their proximity to spliceacceptor/splice donor (SA/SD) regions, whether they were upstream of (orwithin) annotated functional domains, whether they were upstream of (oradjacent to) documented human mutations, as well as their predictedontarget cutting efficiency score (the maximum score is 1.0; higherscore=greater activity) and off target cutting score (the maximum scoreis 100.0; higher score=greater specificity) (Doench et al., NatureBiotechnology, 2014; Hsu et al., 2013). b, Summary of data frompublished studies in which these tumor suppressor genes were inactivatedin the context of Kras^(G12D)-driven lung cancer models c, Each vectorhas a unique sgID and was diversified with random barcodes. The sgID foreach of the vectors and the estimated number of barcodes associated witheach sgRNA is indicated. d, Schematic of the experiment to assess theinitial representation of each sgRNA within Lenti-sg TS-Pool/Cre. e, Thepercent of each sgRNA within Lenti-sg TS-Pool/Cre, as determined bysequencing of samples from three replicate infections. Mean+/−SD isshown. The percent of each vector in the pool deviated only slightlyfrom the expected representation of each vector (red dashed line).

FIG. 13. In vitro sgRNA cutting efficiency. a, Schematic of theexperiment to assess the in vitro cutting efficiency of each sgRNA byinfecting Cas9 cells with lentivirus carrying each individual sgRNA. Wetested three individual sgRNAs for each targeted loci and we report thecutting efficiency of the best sgRNA. b, Cutting efficiency of the bestsgRNA for each targeted tumor suppressor. Cutting efficiency wasassessed by Sanger sequencing and TIDE analysis software (Brinkman etal., Nucl. Acids Res., 2014). c, Schematic of the experiment to assessthe in vitro cutting efficiency of each sgRNA by infecting Cas9 cellswith Lenti-sg TS-Pool/Cre. Cells were harvested 48 hours after infection(transduction), genomic DNA was extracted, the 14 targeted regions werePCR amplified, and the products were sequenced. By calculating the % ofindels at each region, and normalizing to both the representation in thepool and Setd2 indel %, a relative cutting efficiency was determined foreach sgRNA within the pool. d, Relative cutting efficiency of each sgRNAincluding the inert Neo-targeting controls.

FIG. 14. Identification and validation of tumor suppressors at multipletime points using Tuba-seq. a, Percent representation of eachLenti-sgRNA/Cre vector in KT mice 12 weeks after tumor initiation(calculated as 100 times the number of reads with each sgID/all sgIDreads). As there is no Cas9-mediated gene inactivation in KT mice, thepercent of each sgID in these mice represents the percent of viralvectors with each sgRNA in the Lenti-sg TS-Pool/Cre pool. b, Analysis ofthe relative tumor sizes in KT mice (which lack Cas9) 12 weeks aftertumor initiation with Lenti-sg TS-Pool/Cre identified essentiallyuniform tumor size distributions. Relative tumor size at the indicatedpercentiles represents merged data from 10 mice, normalized to theaverage of sg Inert tumors. 95% confidence intervals are shown.Percentiles that are significantly different from sg Inert are in color.c, Estimates of mean tumor size, assuming a lognormal tumor sizedistribution, showed expected minor variability in KT mice.Bonferroni-corrected, bootstrapped p-values are shown. p-values <0.05and their corresponding means are bold. d, Percent representation ofeach Lenti-sgRNA/Cre vector in KT; Cas9 mice 12 weeks after tumorinitiation (calculated as 100 times the number of reads with eachsgID/all sgID reads). e, Tumor sizes at the indicated percentiles foreach sgRNA relative to the average of sg Inert-containing tumors at thesame percentiles. Merged data from 3 KT; Cas9 mice 15 weeks after tumorinitiation with Lenti-sg TS-Pool/Cre is shown. Dotted line represents nochange from Inert. Error bars represent 95% confidence intervals.Percentiles in which the confidence intervals do not overlap the dottedline are in color. f, Estimates of mean tumor size, assuminglognormality, identified sgRNAs with significant growth advantage in KT;Cas9 mice. Bonferroni-corrected, bootstrapped p-values are shown.p-values <0.05 and their corresponding mean estimates are in bold.

FIG. 15. Identification of p53-mediated tumor suppression in KT; Cas9mice with Lenti-sgTS/Cre initiated tumors at two independent timepoints. a,b, Analysis of the relative tumor sizes in KT; Cas9 mice 12weeks (a) and 15 weeks (b) after tumor initiation with Lenti-sgTS-Pool/Cre identify p53 as a tumor suppressor using power-lawstatistics at both time points. Relative tumor size at the indicatedpercentiles is merged data from 8 and 3 mice, respectively, normalizedto the average of sg Inert tumors. 95% confidence intervals are shown.Percentiles that are significantly larger from sg Inert are in color.Power-law p-values are indicated. Note that in this experimental settingonly the very largest sgp53 initiated tumors are greater in size thanthe sg Inert tumors. This is likely partially explained by therelatively poor cutting efficiency of sgp53 (FIG. 13d ). c-f, Percent ofeach size indel at the p53 locus (from ten nucleotide deletions (−10) tothree nucleotide insertions (+3)) were calculated by dividing the numberof reads with indels of a given size by the total number of reads withindels. Inframe indels are shown in grey. We assessed the spectrum ofindels at the p53 locus generated in vitro, in a Cas9 expressing cellline infected (transduced) with Lenti-sg TS-Pool/Cre 48 hours afterinfection (transduction). (c) There is no preference for out of framemutations. We then analyzed three individual KT; Cas9 mice with Lenti-sgTS-Pool/Cre initiated tumors after 15 weeks of disease progression(d-f). There were fewer in-frame indels (−9, −6, −3 and +3) consistentwith selection for out-of-frame loss-of-function alterations in tumorsthat expand, consistent with the tumor suppressive function of p53.These types of analyses, while consistent with the Tuba-seq findings,are imprecise relative to the Tuba-seq platform.

FIG. 16. Analysis of tumor size distributions demonstrates that Lkb1 andSetd2 deficiencies are lognormal. a,b, Size of tumors at the indicatedpercentile (% ile) with sg Lkb1 (a) or sg Setd2 (b) versus sgInert-initiated tumor size at the same percentile. Each percentile wascalculated using all tumors with each sgRNA from all KT; Cas9 mice withLenti-sg TS-Pool/Cre initiated tumors analyzed 12 weeks after tumorinitiation (N=8 mice). The size relative to sg Inert-initiated tumors isindicated with dashed lines. c, Probability density plot for tumorsinitiated with Lenti-sg Setd2/Cre in KT; Cas9 mice with Lenti-sgTS-Pool/Cre initiated tumors shows lognormally distributed tumor sizesvery similar to those seen in KLT mice. This indicates that Setd2deficiency drives tumor growth without providing a significant increasein the generation of, or tolerance to, additional advantageousalterations.

FIG. 17. Confirmation of on-target sgRNA effects. a,b, Percent of eachindel (from ten nucleotide deletions (−10) to four nucleotide insertions(+4)) were calculated by dividing the number of reads with indels of agiven size by the total number of reads with indels within each toptumor suppression gene. (a) Average percentage and standard deviation ofthree KT; Cas9 mice with Lenti-sg TS-Pool/Cre-initiated tumors are shownfor Setd2, Lkb1, Rb1, and the average of the three targeted sites in Neo(Neo1-3). Inframe mutations are shown in grey. Average and standarddeviations for Neo1-3 was calculated by averaging all three mice and allthree Neo target sites as a single group. In general, there were fewerinframe indels (−9, −6, −3 and +3) consistent with selection forout-of-frame loss-of-function alterations in these genes in tumors thatexpand. (b) We also assessed the spectrum of indels generated in vitro,in a Cas9-expressing cell line infected (transduced) with Lenti-sgTS-Pool/Cre 48 hours after infection (transduction). We detected nopreference for inframe mutations in any of these genomic locations,suggesting that the bias in the KT; Cas9 mice is most likely due toadvantageous expansion of tumors with out-of-frame indels (i.e., nullallele). c, Kaplan-Meier survival curve of KT and KT; Cas9 mice withLenti-sg Smad4/Cre-induced tumors. CRISPR/Cas9-mediated inactivation ofSmad4 in the presence of oncogenic Kras^(G12D) does not reduce survival,suggesting limited, if any, increase in tumor growth from Smad4inactivation. d, The majority of tumors in Lenti-sg Smad4/Cre infected(transduced) KT; Cas9 mice had lost Smad4 protein expression compared toKT mice infected (transduced) with the same virus, consistent with indelcreation at the Smad4 locus. Scale bars=50 μm. e, Several tumors inLenti-sg TS-Pool/Cre-infected (transduced) KT; Cas9 mice had a distinctpapillary histology, uniformly large nuclei, and were Sox9 positive,consistent with the published phenotype of Apc-deficient, Kras-drivenlung tumors (Sanchez-Rivera et al., Nature, 2014). RepresentativeSox9-negative and Sox9-positive tumors are shown. Scale bars=100 μm(top) and 25 μm (bottom).

FIG. 18. Additional images showing increased tumor burden in mice withCRISPR/Cas9-mediated inactivation of Setd2 using each of two independentsgRNAs. Additional representative fluorescence dissecting scope imagesof lung lobes from KT; Cas9 mice with tumors initiated withLenti-sgNeo2/Cre (left), Lenti-sg Setd2#1/Cre (middle), or Lenti-sgSetd2#2/Cre (right) analyzed 9 weeks after tumor initiation. Lung lobesare outlined with white dashed lines. Scale bars=5 mm.

FIG. 19. Comparison of systems to assess tumor suppressor gene functionin lung adenocarcinoma mouse models. The method of tumor suppressor geneinactivation (Cre/LoxP-mediated deletion of a floxed allele versusCRISPR/Cas9-mediated genome editing), the ability to quantify tumornumber and size through genetic barcoding of individual tumors, and theability to inactivate multiple genes in a pooled format is indicated.Particularly relevant advantages and disadvantages of each system areshown, as well as example references. All highlighted studies are inlung cancer except Maresch et al. who used pooled sgRNA transfection tostudy pancreatic cancer. The reality of using floxed alleles to assesstumor suppressor gene function in lung adenocarcinoma models is bestexemplified by the fact that over the past 15 years only six of thetumor suppressor genes that we queried have been investigated usingfloxed alleles in combination with Kras^(LSL-G12D). The lack ofquantitative methods also severely hampers the identification of geneswith only moderate tumor suppressive effects due to known and unknowntechnical and biological variables (e.g. reproducibility of tumorinitiation, gender, age, and strain of mice). Data generated by deletinggenes with floxed alleles is also limited by the difficulty in comparingbetween different experimental setups used in different laboratories(e.g. different viral titer, time after initiation, method ofquantification, mouse strain). Thus the relative effect of differenttumor suppressor genes is difficult to glean from the literature.Finally, the quantification of individual tumor cell number by tumorbarcoding provides not only unprecedented precision but also uncoversgene-specific effects on tumor size distributions that likely reflectdistinct functional mechanisms of tumor suppression.

FIG. 20. Statistical properties of tumor size distributions and thecovariance of sgRNA tumor sizes across mice. a. The mean and variance ofeach sgID distribution in every mouse with Lenti-sg Pool/Cre initiatedtumors. Mouse genotypes are colored as indicated. In general, varianceincreased with the square of the mean for all genotypes, suggesting thata log-transformation of lesion size should stabilize variance and avoidheteroskedasticity. Some distributions exhibit a variance that increasedby more than the square of the mean. b-d. Mouse-to-mouse variability inresponse to genetic alterations was interrogated in KT; Cas9 micesacrificed at 12 weeks. The covariance of the LN MLE mean of each sgRNAin each mouse was investigated. Genotype means sizes positivelycorrelated with each other across mice (e.g. a mouse with larger sg Lkb1tumors also harbored larger sg Setd2 tumors.) PCA decomposition of thecorrelation matrix amongst all 12 sgRNAs (sg Inerts consolidated)uncovered a substantial level of mouse-to-mouse variability explicableby a single Principle Component (PC1) vector. Each dot represents asingle mouse projected onto PC1, which explains 75% of observedvariability between mice in sgRNA mean sizes. (b) PC1 correlates withoverall lung weight and (c) mean lesion size, indicating that mice withlarger tumors are more susceptible to tumor growth driven by strongdrivers (PC1 correlated with sg Setd2 and sg Lkb1 size, data not shown.)(d) The mice do not appear to form distinct clusters when projected ontothe first two Principle Components. Replicate mice were almost alwayssiblings housed in the same cages. We minimized extrinsic sources ofnoise using a Mixture of Principal Components model (see Methods.)

FIG. 21. Mathematical models of tumor progression.

FIG. 22. Frequency of lentiviral infections (transductions) compared tosize difference between each lesion and its nearest neighbor in the samemouse.

FIG. 23. A platform that integrates AAV/Cas9-mediated somatic HDR withtumor barcoding and sequencing to enable the rapid introduction andfunctional investigation of putative oncogenic point mutations in vivo.a-d. Schematic overview of the pipeline to quantitatively measure the invivo oncogenicity of a panel of defined point mutations. A library ofAAV vectors was generated such that each AAV contains 1) a template forhomology directed repair (HDR) containing a putatively oncogenic pointmutation and a random DNA barcode encoded in the adjacent wobble bases,2) an sgRNA targeting the endogenous locus for HDR, and 3)Cre-recombinase to activate a conditional Cas9 allele(H11^(LSL-Cas9) and other Cre-dependent alleles in genetically engineered mice (a). The AAV library is delivered to a tissue of interest (b). Following transduction, a subset of cells undergo AAV/Cas)9-mediatedHDR in which the locus of interest is cleaved by Cas9 at the sgRNAtarget site and repaired using the AAV HDR template. This results in theprecise introduction of the desired point mutation and a unique DNAbarcode into the targeted locus (c). Somatic cells engineered with apoint mutation may develop into de novo tumors if the introducedmutation is sufficient to initiate tumorigenesis and drive tumor growth.d, Two independent approaches can be used to analyze tumors: 1) tumorscan be sequenced individually to characterize both alleles of thetargeted gene, or 2) barcoded mutant HDR alleles from entire bulktumor-bearing tissues can be deep sequenced to quantify the number andsize of tumors with each mutation. e. AAV vector pool for Cas9-mediatedHDR into the endogenous Kras locus (AAV-Kras^(HDR)/sg Kras/Cre). Eachvector contains an HDR template with 1 of 12 non-synonymous Krasmutations at codons 12 and 13 (or wild type Kras), silent mutationswithin the PAM and sgRNA homology region (PAM*), and an 8-nucleotiderandom barcode within the wobble positions of the downstream codons forDNA barcoding of individual tumors. f. Representation of each Kras codon12 and 13 allele in the AAV-Kras^(HDR)/sg Kras/Cre plasmid library. g.Diversity of the barcode region in the AAV-Kras^(HDR)/sg Kras/Creplasmid library.

FIG. 24. AAV/Cas9-mediated somatic HDR initiates oncogenic Kras-drivenlung tumors that can progress into a metastatic state. a. Schematic ofthe experiment to introduce point mutations and a DNA barcode into theendogenous Kras locus of lung epithelial cells in Rosa26^(LSL-tdTomato);H11^(LSL-Cas9) (T; H11^(LSL-Cas9)), p53^(flox/flox); T; H11^(LSL-Cas9)(PT; H11^(LSL-Cas9)), and Lkb1^(flox/flox); T; H11^(LSL-Cas9) (LT;H11^(LSL-Cas9)) mice by intratracheal administration ofAAV-Kras^(HDR)/sg Kras/Cre. b. Representative images ofTomato^(positive) lung tumors and histology in AAV-Kras^(HDR)/sgKras/Cre-treated LT; H11^(LSL-Cas9); PT; H11^(LSL-Cas9), and T;H11^(LSL-Cas9) mice. Scale bars=5 mm. c. Quantification of lung tumorsin the indicated genotypes of mice infected (transduced) with theindicated AAV vectors (with and without sg Kras). Each dot representsone mouse. Kras^(LSL-G12D); LT (KLT) and Kras^(LSL-G12D); PT (KPT) micetransduced with a 1:10,000 dilution of AAV-Kras^(HDR)/sg Kras/Credeveloped approximately half as many tumors as the PT; H11^(LSL-Cas9)and LT; H11^(LSL-Cas9) mice infected (transduced) with undiluted virus.Thus, assuming that all Kras^(HDR) alleles in the AAV-Kras^(HDR)/sgKras/Cre library are oncogenic, this suggests that AAV/Cas9-mediated HDRoccurs in approximately 0.02% of transduced cells. Alternatively, ifonly 20% of the mutant alleles in the AAV-Kras^(HDR)/sg Kras/Cre libraryare assumed to drive tumor formation, then the rate of HDR isapproximately 0.1%. d. Representative FACS plot showingTomato^(positive) disseminated tumor cells (DTCs) in the pleural cavityof an LT; H11^(LSL-Cas9) mouse with AAV-Kras^(HDR)/sg Kras/Cre-initiatedlung tumors. e. Histology of a metastasis from an AAV-Kras^(HDR)/sgKras/Cre-initiated lung tumor in a PT; H11^(LSL-Cas9) mouse. Scalebar=50 μm. f. Diverse HDR-generated oncogenic Kras alleles in individuallung tumors. Number of tumors with each allele is indicated. Allelesthat were not identified in any lung tumors are not shown.

FIG. 25. Introduction of mutant Kras variants into somatic pancreas andmuscle cells by AAV/Cas9-mediated HDR drives the formation of invasivecancers. a. Schematic of retrograde pancreatic ductal injection ofAAV-Kras^(HDR)/sg Kras/Cre into PT; H11^(LSL-Cas9) mice to inducepancreatic cancer. b. Histology of pancreatic tumors initiated byretrograde pancreatic ductal injection of AAV-Kras^(HDR)/sg Kras/Creinto PT; H11^(LSL-Cas9) mice. Scale bars=75 μm. c. Histology ofmetastases in the lymph node (upper panel) and diaphragm (lower panel)in PT; H11^(LSL-Cas9) mice with primary PDAC. Scale bars=50 μm. d.HDR-generated oncogenic Kras alleles in pancreatic tumor masses. Numberof tumors with each allele is indicated. Alleles that were notidentified in any pancreatic tumor masses are not shown. e. Schematic ofintramuscular injection of AAV-Kras^(HDR)/sg Kras/Cre into thegastrocnemii of PT; H11^(LSL-Cas9) mice to induce sarcomas. f,g.Histology of stereotypical sarcoma (f) and invasive sarcoma (g)initiated by intramuscular injection of AAV-Kras^(HDR)/sg Kras/Cre intothe gastrocnemii of PT; H11^(LSL-Cas9) mice. Scale bars=75 μm. h.HDR-generated oncogenic Kras alleles in sarcomas. Number of tumors witheach allele is indicated. Alleles that were not identified in anysarcomas are not shown. These data document clonal marking of celllineages across multiple tissues.

FIG. 26. Multiplexed, quantitative analysis of Kras mutant oncogenicityusing AAV/Cas9-mediated somatic HDR and high-throughput sequencing ofindividually barcoded tumors. a. Pipeline to quantitatively measureindividual tumor size and number from bulk lung samples byhigh-throughput sequencing of tumor barcodes. b. Number of lung tumorsharboring each mutant Kras allele normalized to its initialrepresentation (mutant representation in the AAV plasmid library/WTrepresentation in the AAV plasmid library) and relative to WT (mutanttumor #/WT tumor #). Variants present in significantly more tumors thanWT (p<0.01) are colored blue; darker blue indicates no significantdifference from G12D (p>0.05), lighter blue indicates significantly lesstumors with that variant than G12D (p<0.01). c. p-values from atwo-sided multinomial chi-squared test of the number of lung tumors witheach Kras variant across different genotypes. Significant p-values(p<0.05) are bold. d,e. Lung tumor size distributions for Kras variantsidentified as oncogenic in b across all LT; H11^(LSL-Cas9) (d) or PT;H11^(LSL-Cas9) (e) mice. Each dot represents one tumor with a uniqueKras variant-barcode pair. The size of each dot is proportional to thesize of the tumor it represents, which is estimated by normalizing tumorread counts to the normalization control reads counts. f. DiverseHDR-generated Kras alleles identified by tumor barcode sequencing ofpancreatic tumor masses. Number of uniquely barcoded tumors with eachallele is indicated. Alleles that were not identified in any pancreastumor masses are not shown. g. High-throughput sequencing of the primarypancreatic tumor mass and metastases from a single AAV-Kras^(HDR)/sgKras/Cre-treated PT; H11^(LSL-Cas9) mouse uncovered a diverse spectrumof mutant Kras alleles and enabled the establishment of clonalrelationships between primary tumors and their metastatic offspring.Each dot represents one tumor with the indicated Kras variant and aunique barcode within the indicated sample. Dots that are linked by acolored line harbor the same barcode, suggesting that they are clonallyrelated. The size of each dot is scaled according to the size of thetumor it represents (diameter of the dot=relative size^(1/4)). Since thesize of pancreatic tumors is not normalized to a control, tumors sizescan only be compared within the same sample. Thus, the largest tumor ineach sample is set to the same standard size.

FIG. 27. Design, generation, and validation of an AAV library formultiplexed mutation of Kras. a. Sequence of the three sgRNAs targetingKras exon 2. Cutting efficiency of each sgRNA was determined bysequencing DNA from Cas9-expressing MEFs 48 hours after transductionwith lentiviral vectors encoding each sgRNA. All three sgRNAs inducedindel formation at the targeted loci. Thus, the sgRNA targeting thesequence closest to Kras codons 12 and 13 (sg Kras#3) was used for allsubsequent experiments to increase the likelihood of HDR. b. Synthesizedlibrary of dsDNA fragments containing wild type (WT) Kras sequence pluseach of the 12 non-synonymous, single nucleotide Kras mutants at codons12 and 13, silent mutations within the PAM and sgRNA homology region(PAM*), and an 8-nucleotide random barcode within the wobble positionsof the downstream codons for barcoding of individual tumors. Each Krasallele can be associated with ˜2.4×10⁴ unique barcodes. Fragments alsocontained restriction sites for cloning. c. AAV vector library wasgenerated by massively ligating synthesized regions into a parental AAVvector creating a barcoded pool with WT Kras and all 12single-nucleotide, non-synonymous mutations in Kras codons 12 and 13. d.Position of Kras exon 2 within the Kras^(HDR) template. The lengths ofthe homology arms are shown. e. Schematic of the experiment to test forHDR bias. A Cas9-expressing cell line was transduced withAAV-Kras^(HDR)/sg Kras/Cre and then sequenced to quantify HDR events. f.Schematic of the PCR strategy to specifically amplify Kras^(HDR) allelesintroduced into the genome via HDR. Forward primer 1 (F1) binds to thesequence containing the 3 PAM* mutations, while reverse primer 1 (R1)binds the endogenous Kras locus, outside the sequence present in thehomology arm of the Kras^(HDR) template. F2 binds to the IIluminaadaptor added by F1, R2 binds to a region near exon 2, and R3 binds tothe Illumina adapter added in the same reaction by R2. g. Representationof each Kras allele within the endogenous Kras locus generated throughHDR in Cas9-expressing cells in culture transduced with theAAV-Kras^(HDR)/sg Kras/Cre vector library. h. Frequency of HDR eventsfor each Kras^(HDR) allele plotted against the initial frequency of eachKras mutant allele in the AAV-Kras^(HDR)/sg Kras/Cre plasmid libraryused to generate the viral library. High-correlation between the initialplasmid library and the representation of mutant Kras alleles followingHDR suggests little to no HDR bias.

FIG. 28. Identification of an optimal AAV serotype for adult lungepithelial cell transduction. a. Outline of the experiment to screen 11AAV serotypes for adult lung epithelial cell transduction. An AAV vectorencoding GFP was packaged with different AAV capsid serotypes andadministered intratracheally to wild-type recipient mice. 5 dayspost-treatment, the lungs were dissociated and the percent ofGFP^(positive) epithelial cells was determined by flow cytometry. b.Different AAV serotypes can be produced at different concentrations. Ourgoal was to identify the AAV serotypes capable of delivering DNAtemplates to lung epithelial cells, which is largely dictated by boththe achievable viral titer and the per virion transduction efficiency.Thus, we did not normalize the titer of the AAV serotypes beforeinfection (transduction), but rather determined the percent infection(transduction) following administrations of 60 μl of undiluted, purifiedvirus. c. To assess the percent of lung epithelial cell transduced bythe different AAV serotypes, we dissociated lungs of infection(transduction) mice into single cell suspensions and performed flowcytometry for GFP as well as for markers of hematopoietic cells (CD45,Ter119, and F4/80), endothelial cells (CD31), and epithelial cells(EpCAM). Plots show FSC/SSC-gated, viable (DAPI^(negative)), lungepithelial (CD45/Ter119/F4-80/CD31negative, EpCAM^(positive)) cells. Thepercent GFP^(positive) epithelial cells in each sample is indicatedabove the gate. AAV8, AAV9, and AAVDJ were considerably better than allother serotypes (including AAV6 which failed to lead to efficient HDR inPlatt et al., Cell, 2014), consistent with the high maximal titers ofthese serotypes. We chose to use AAV8 based on this data and thedocumented ability of AAV8 to efficiently transduce many other mousecell types in vivo.

FIG. 29. AAV/Cas9-mediated in vivo HDR in lung epithelial cellsinitiates primary tumors that can progress to gain metastatic ability.a. Schematic of the experiment to introduce point mutations into theendogenous Kras locus and barcode lung epithelial cells inLkb1^(flox/flox), R26^(LSL-Tomato); H11^(LSL-Cas9) (LT; H11^(LSL-Cas9)),p53^(flox/flox); R26^(LSL-Tomato); H11^(LSL-Cas9) (PT; H11^(LSL-Cas9))and R26^(LSL-Tomato); H11^(LSL-Cas9) (T; H11^(LSL-Cas9)) mice byintratracheal administration of AAV-Kras^(HDR)/sg Kras/Cre. b. Lightimages that correspond to the fluorescence images in FIG. 2a . Highermagnification histology images document adenocarcinoma histology andgreater nuclear atypia in the p53-deficient tumors. Upper scale bars=5mm. Lower scale bars=50 μm. c. Additional examples of AAV-Kras^(HDR)/sgKras/Cre-induced lung tumors in LT; H11^(LSL-Cas9), PT; H11^(LSL-Cas9),and T; H11^(LSL-Cas9) mice. Scale bars=5 mm. Note that, due to the hightransduction efficiency, most lung cells express Tomato, but the tumorsare much brighter because of the large number and density of cells ineach tumor. d. Total lung weight in mice of each genotype with tumorsinitiated with AAV-Kras^(HDR)/sg Kras/Cre. Each dot represents onemouse. e. Number of surface lung tumors identified under a fluorescencedissecting scope in mice of each genotype infected (transduced) withAAV-Kras^(HDR)/sg Kras/Cre diluted 1:10. Each dot represents one mouse.f. Histology of a lymphatic micrometastasis that formed in a PT;H11^(LSL-Cas9) mouse with AAV-Kras^(HDR)/sg Kras/Cre-initiated lungtumors. Scale bar=50 μm. g. Number of mice of each genotype that haddisseminated tumors cells in the pleural cavity (DTCs>10) and lymph nodemetastases. The numbers represent the number of mice with DTCs ormetastases/total number of mice analyzed.

FIG. 30. Nuclease-free AAV-mediated HDR does not occur at a high enoughrate to initiate large numbers of lung tumors. a. Schematic of controlAAV vector library that contains a 2.5 kb Kras HDR template with the 12single-nucleotide, non-synonymous mutations and barcode, but without thesgRNA targeting Kras. b. Representation of each Kras codon 12 and 13allele in the AAV-Kras^(HDR)/Cre plasmid pool. Percentages are theaverage of triplicate sequencing. c. Titer of the AAV vector libraries(vg=vector genomes). Importantly, the control AAV-Kras^(HDR)/Cre viralpreparation is higher titer than AAV-Kras^(HDR)/sg Kras/Cre. d.Quantification of the number of LT, PT, and T mice that developed tumorsafter administration of 60 μL of undiluted or 1:10 dilutedAAV-Kras^(HDR)/Cre pool.

FIG. 31. Analysis of individual tumors identifies oncogenic Kras allelesand uncovers indels in the non-HDR Kras allele. a. Example sequencingtrace of a Kras^(HDR) allele with PAM* mutations, a G12D mutation, and abarcode. b. Sequences of four representative oncogenic Kras allelesdetected in individual lung tumors by Sanger sequencing. Each primarytumor analyzed had a unique variant-barcode pair, as expected given˜2.4×10⁴ possible barcodes per variant. The altered bases in theAAV-Kras^(HDR) template sequence and the wild type Kras sequence at thislocus are shown for reference. c. HDR events generally occurred outsideof the two engineered restriction sites. However, some tumors had Krasalleles consistent with recombination between exon 2 and one of therestriction sites, suggesting recombination very close to the Cas9/sgKras-induced double-strand DNA break. d. Diagram of oncogenic Krasalleles in individual tumors that did not undergo perfect HDR. Bothperfect and imperfect HDR events are found in each mouse genotype(perfect HDR in 14/30 tumors in LT; H11^(LSL-Cas9) mice and 3/7 tumorsin PT; H11^(LSL-Cas9) mice). Imperfect HDR events included alleleslikely integrating into the Kras locus through homologous recombinationof the 5′ end of the AAV-Kras^(HDR) template upstream of exon 2 andligation of the 3′ end of the AAV-Kras^(HDR) template to the exon 2region immediately downstream of the Cas9/sg Kras-induced double-strandDNA break. This imperfect HDR resulted in insertions or deletions in theintronic sequence downstream of Kras exon 2. Insertions and deletionswere variable in length (sizes approximated by Sanger sequencing or gelelectrophoresis) and sometimes included part or all of the wild typeexon 2, or in rare cases, segments of the AAV-Kras^(HDR)/sg Kras/Crevector. None of these partial HDR events were predicted to altersplicing from the mutant exon 2 to exon 3, consistent with therequirement for expression of the oncogenic Kras allele for tumorformation. e,f. The oncogenic Kras allele in large individual tumorsfrom treated PT; H1^(LSL-Cas9) and LT; H11^(LSL-Cas9) mice was almostalways accompanied by inactivation of the other Kras allele throughCas9-mediated indel formation in exon 2. Sanger sequencing identifiedindels adjacent to the PAM sequence in 47/48 (98%) of individual tumors.Example indels (e) and a summary of all indels (f) are shown. NDindicates that a wild type allele could not be detected, which isconsistent with either loss of heterozygosity, a very large indel, or alarge deletion that encompassed one of the primer binding sites.

FIG. 32. HDR-mediated introduction of oncogenic mutations into theendogenous Kras locus in pancreatic cells leads to the formation ofpancreatic ductal adenocarcinoma. a. Schematic of retrograde pancreaticductal injection of AAV-Kras^(HDR)/sg Kras/Cre into PT; H11^(LSL-Cas9)mice to induce pancreatic cancer. b. Representative light andfluorescence images of pancreatic tumors that developed in PT;H11^(LSL-Cas9) mice transduced with AAV-Kras^(HDR)/sg Kras/Cre. Scalebars=5 mm. c. Histology images of different stages of pancreatic tumorprogression including a pre-cancerous PanIN lesion (upper left), awell-differentiated tumor region (top right), and poorly differentiatedPDAC (bottom left). Bottom right shows the development of acollagen-rich stromal environment (stained with Trichrome) within PDAC.Scale bars=75 μm. d. Representative FACS plots showing Tomato^(positive)disseminated tumor cells (DTCs) in the peritoneal cavity of a PT;H11^(LSL-Cas9) mouse with AAV-Kras^(HDR)/sg Kras/Cre-initiated PDAC.Plot shows FSC/SSC-gated viable cancer cells(DAPI/CD45/CD31/F4-80/Ter119^(negative)). e. HDR-induced PDACs canprogress to gain metastatic ability, seeding metastases in lymph nodesand on the diaphragm. Light and fluorescence dissecting scope images areshown. Scale bars=3 mm. f. Incidence of PDAC, DTCs in the peritonealcavity, and metastases in the indicated genotypes of mice (shown as thenumber of mice with cancer, DTCs, or metastases out of the total numberof mice analyzed) 3-13 months post-infection (transduction) with theindicated AAV vector libraries.

FIG. 33. HDR-mediated induction of oncogenic Kras in skeletal muscleinduces sarcomas. a. Schematic of intramuscular injection ofAAV-Kras^(HDR)/sg Kras/Cre into the gastrocnemii of PT; H11^(LSL-Cas9)mice to induce sarcomas. b. Representative whole mount light (top panel)and fluorescence dissecting scope (bottom panel) images of mousegastrocnemii following injection with AAV-Kras^(HDR)/sg Kras/Cre. Rightgastrocnemius has sarcoma, while the left does not, despite efficienttransduction as evidenced by widespread Tomatopositive tissue (data notshown). Scale bars=5 mm. c. Images of histological H&E sectionsconfirming the presence of sarcoma with stereotypical histology and alsoinvasion into the surrounding muscle. Scale bars=75 μm. d. Incidence ofsarcomas in PT; H11^(LSL-Cas9) mice 3-7 months after intramuscularinjection of AAV-Kras^(HDR)/sg Kras/Cre. Incidence represents the numberof mice that developed sarcomas out of the total number of miceinjected. One of the 7 treated mice has not yet been analyzed but didnot have an obvious sarcoma six months post-infection (transduction). e.Sequencing of the Kras^(HDR) locus in a sarcoma reveals a mutant Krasallele and barcode.

FIG. 34. Samples and preparation for Illumina® sequencing of bulk lungtissue to quantify the size and number of lung tumors with each mutantKras allele. a. Bulk lung tissue samples from mice intratracheallyadministered with AAV-Kras^(HDR)/sg Kras/Cre for Illumina® sequencing ofbarcoded Kras^(HDR) alleles. Sample name, mouse genotype, and dilutionof AAV-Kras^(HDR)/sg Kras/Cre are indicated. The weight, tumor number,number of dissected tumors, as well as the amount of DNA amplified andthe number of PCR reactions pooled for Illumina® sequencing for eachsample are shown. Repeat samples are technical replicates. ND=No data.b. Simplified pipeline for the normalization of sequencing reads frombulk lung samples using reads from a benchmark control of known cellnumber to enable estimation of cell number in each tumor and allow datafrom separate mice to be combined.

FIG. 35. Reproducibility of barcode sequencing-based parallel analysisof tumor genotype, size, and number from bulk tissue. a-d. Regressionplot of individual tumors with the indicated Kras^(HDR) allele and aunique barcode detected by high-throughput sequencing across technicalreplicates (i.e. independent DNA extraction from bulk tissue lysate andPCR reactions). Replicates in a and b were PCR amplified using primerswith different multiplexing tags, but were run on the same sequencinglane. Replicates in c and d were PCR amplified using the same primers,but were run on different sequencing lanes. Mice with above averagetumor burden (a,c) and below average tumor burden (b,d), as estimatedmeasured by bulk lung weight, were analyzed to confirm the technical andcomputational reproducibility of this pipeline across samples ofvariable tumor number.

FIG. 36. High-throughput barcode sequencing of tumors from bulk lungtissue uncovers diverse numbers and sizes of tumors. a-c. Tumor sizedistributions of all Kras variants across all LT; H11^(LSL-Cas9) (N=6)(a), PT; H11^(LSL-Cas9) (N=7) (b) or T; H11^(LSL-Cas9) (N=3) (c) mice.Each dot represents a tumor with a unique Kras variant-barcode pair. Thesize of each dot is proportional to the size of the tumor it represents,which is estimated by normalizing tumor read counts to the normalizationcontrol reads counts. Lesions harboring WT Kras^(HDR) alleles arethought to be hitchhikers in tumors with oncogenic Kras^(HDR) alleles(see Methods). d,e. Tables of raw (d) and normalized (e) number oftumors harboring each Kras variant across each genotype (includingtumors with each variant that were identified by individual tumordissection and analysis). In e, the number of tumors harboring each Krasvariant is normalized to the initial representation of each variant inthe AAV plasmid library and to the number of lesions harboring a WTallele within the same genotype. Note that the color intensity scale ofthe heatmaps in e is unique to each genotype for ease of comparison.

FIG. 37. High-throughput sequencing of pancreatic tumor masses andmetastases identifies oncogenic Kras mutants. a. Bulk pancreas tissueand metastasis samples from mice administered with AAV-Kras^(HDR)/sgKras/Cre by retrograde pancreatic ductal injection for Illuminasequencing of barcoded Kras^(HDR) alleles. Sample name, mouse genotype,viral dilution, and tissue are indicated. The Kras^(HDR) alleles presentin distinct regions of the primary tumor masses as well as metastaseswere analyzed by Illumina® sequencing after FACS isolating FSC/SSC-gatedviable cancer cells (DAPI/CD45/CD31/F4-80/Ter119^(negative)) from thesesamples. b. Analysis pipeline to identify Kras^(HDR) alleles inAAV-Kras^(HDR)/sg Kras/Cre-initiated tumor masses within the pancreataof PT; H11^(LSL-Cas9) mice. c. Multi-region sequencing of a largepancreatic tumor mass in a single AAV-Kras^(HDR)/sg Kras/Cre-treated PT;H11^(LSL-Cas9) mouse uncovered a diverse spectrum of mutant Kras allelesand linked primary tumors with their metastatic offspring. Each dotrepresents a tumor with the indicated Kras variant and a barcode uniqueto the indicated sample (labeled 1-4). Dots connected across differentprimary tumor samples (labeled 1-3) shared the same Kras variant-barcodepair, and are thus presumably regions of the same primary tumor thatwere present in multiple samples. A colored line link primary tumors andlymph node metastases harboring the same Kras variant-barcode pair,indicating a clonal relationship. The size of each dot is scaledaccording to the size of the tumor it represents (diameter of thedot=relative size^(1/2)). Since the size of pancreatic tumors is notnormalized to a control, tumor sizes can only be compared to othertumors within the same sample. Thus, the largest tumors within eachsample have been scaled to the same standard size. g=gallbladder,sto=stomach, duo=duodenum, pan=pancreas, sp=spleen, ln=mesenteric lymphnodes.

FIG. 38. Relationship between the in vivo oncogenicities and biochemicalbehaviors of Kras mutants. a-c. Relative number of lung tumors in micetransduced with AAV-Kras^(HDR)/sg Kras/Cre (see FIG. 4b ) as a functionof the indicated biochemical property reported in Hunter et al., 2015.Relative lung tumor number is normalized to the initial representationof each Kras variant in the AAV-Kras^(HDR)/sg Kras/Cre plasmid pool.Vertical bars represent the 95% confidence interval for the normalizedrelative lung tumor number. Horizontal bars represent the standard errorof the mean of three replicate experiments as described in Hunter etal., 2015. P120GAP was used to determine GAP-stimulated GTP hydrolysisrates (Hunter et al., 2015). d-f. Number of pancreatic tumors in micetransduced with AAV-Kras^(HDR)/sg Kras/Cre (see FIG. 4f ) as a functionof the indicated biochemical property reported in Hunter et al., 2015.Vertical bars represent the 95% confidence interval for pancreas tumornumber. Horizontal bars represent the standard error of the mean ofthree replicate experiments as described in Hunter et al., 2015. P120GAPwas used to determine GAP-stimulated GTP hydrolysis rates (Hunter etal., 2015).

FIG. 39. Investigating combined genetic alterations: p53 deficiencyalters the growth effects of tumor suppression in KrasG12D-driven lungtumors in vivo. a. Tuba-seq approach to study combinatorial tumorsuppressor inactivation in vivo. Tumors were initiated with Lenti-sgTS-Pool/Cre (containing four inert sgRNA vectors and eleven vectorstargeting known and candidate tumor suppressor genes) in three differentgenetically-engineered mouse backgrounds: Kras^(LSL-G12D/+);Rosa26^(LSL-tdTomato); H11^(LSL-Cas9) (KT; Cas9), KT; p53^(flox/flox);Cas9 (KPT; Cas9), and KT; Lkb1^(flox/flox); Cas9 (KLT; Cas9). Each sgRNAvector contains a unique sgID and a random barcode, which was used toquantify individual tumor sizes via deep sequencing. b. Analysis of therelative tumor sizes in KT; Cas9 mice 15 weeks after tumor initiation.Relative size of tumors at the indicated percentiles is merged data from10 mice, normalized to the average size of sg Inert tumors. Error barsthroughout this study denote 95% confidence intervals determined bybootstrap sampling. Percentiles that are significantly different from sgInert are in color. c. Estimates of mean tumor size, assuming alognormal tumor size distribution, identified sgRNAs that significantlyincreased growth in KT; Cas9 mice. Bonferroni-corrected, bootstrappedP-values are shown. sgRNAs with P-values <0.05 are bold. d,e. Same asb,c, except for merged data from 12 KPT; Cas9 mice. f. Abundance ofindels at targeted loci relative to median of genome-targeting inertsgRNAs Neo1-3. Coloring according to a. g. Functional mutations in TP53and RB1 in human lung adenocarcinomas from TCGA and GENIE datasets(N=1792). RB1 and TP53 alterations co-occur.

FIG. 40. Investigating combined genetic alterations: Attenuated effectsof tumor suppressor inactivation in Lkb1-deficient tumors furtherhighlights a rugged fitness landscape. a. Tumor sizes at the indicatedpercentiles for each sgRNA relative to the average of sgInert-containing tumors at the same percentiles. Merged data from 13 KT;Lkb1^(flox/flox); Cas9 (KLT; Cas9) mice 15 weeks after tumor initiationwith Lenti-sg TS-Pool/Cre is shown. Percentiles that are significantlydifferent from sg Inert are in color. b. Estimates of mean tumor size,assuming a lognormal tumor size distribution, identified sgRNAs thatsignificantly increase growth in KLT; Cas9 mice. Bonferroni-corrected,bootstrapped P-values are shown. sgRNAs with P-values <0.05 are bold. c.Mutual exclusivity of LKB1 (STK11) and SETD2 mutations in human lungadenocarcinomas from TOGA and GENIE datasets (N=1792). d. Tumor sizes inKPT; Cas9 mice with Lenti-sg Setd2/Cre-initiated tumors (N=7) versusKPT; Cas9 mice with Lenti-sgNeo2/Cre initiated tumors (N=3). Lenti-sgSetd2/Cre-initiated tumors have an LN mean that is 2.4 times higher thanLenti-sgNeo2/Cre-initiated tumors and a 95^(th) percentile tumors sizethat is 4.6 times higher. e. Tumor sizes in KLT; Cas9 mice with Lenti-sgSetd2/Cre-initiated tumors (N=7) versus KLT; Cas9 mice withLenti-sgNeo2/Cre initiated tumors (N=5). The relative LN Mean andrelative 95th percentile are 2.2 and 2.8, which are both significantlyless than in FIG. 2d (P<0.04, and P<0.0001 respectively). f. Pearsoncorrelations of fitness effect of tumor suppressors (determined by LNmean) across genetic backgrounds. sgp53 and sg Lkb1 growth rates areexcluded in KPT; Cas9 and KLT; Cas9 mice. *P<0.05, ****P<0.0001. g.Differential effect of each tumor suppressor gene within the context ofoncogenic Kras-driven lung tumors, as well as with coincident p53- orLkb1-deficiency. 95th percentiles that significantly deviate from sgInert tumors are shown in blue. h. Likelihood of identifying candidatetumor suppressors as a driver (as defined in g) versus the number ofgenetic backgrounds studied. All genetic contexts were averaged.

FIG. 41. The current state of genetically-engineered mouse models oflung cancer for the analysis of the putative tumor suppressoralterations in this study and the frequency of these genomic alterationsin human lung adenocarcinoma. a. Summary of data from published studiesin which the putative tumor suppressor genes studied here wereinactivated in the context of oncogenic Kras-driven lung cancer models,with or without inactivation of p53 or Lkb1. b. The percent of tumorswith potentially inactivating alterations (frameshift or non-synonymousmutations, or genomic loss) in each tumor suppressor gene for all tumors(All) as well as for tumors with potentially inactivating alterations inTP53 (TP53^(mut)) or LKB1 (LKB1^(mut)). The percent of tumors with eachtype of alteration is indicated. Data is shown for two clinical cancergenomics studies: The Cancer Genome Atlas (TOGA, 2014), and the GenomicsEvidence Neoplasia Information Exchange (GENIE, 2017) database.

FIG. 42. Description of multiplexed lentiviral vectors, tumorinitiation, and Tuba-seq pipeline to quantify tumor size distributionsin vivo. a. Lenti-sg TS-Pool/Cre contains four vectors with inert sgRNAsand eleven vectors with tumor suppressor gene targeting sgRNAs. EachsgRNA vector contains a unique sgID and a random barcode.NT=Non-Targeting. b. Schematic of the sgID-barcode region of the vectorsin Lenti-sg TS-Pool/Cre. Lenti-sg TS-Pool/Cre contains vectors withfifteen different 8-nucleotide unique identifiers (sgIDs) which link agiven sgID-barcode read to a specific sgRNA. These vectors also containsa 15-nucleotide random barcode element. This double barcode systemallows identification of individual tumors, as well as the sgRNA in thevector that initiates each tumor. c. Transduction of lung epithelialcells with the barcoded Lenti-sg TS-Pool/Cre pool initiates lung tumorsin genetically engineered mouse models with (1) a Cre-regulatedoncogenic KrasG12D (Kras^(LSL-G12D/+)) allele, (2) a Cre reporter allele(Rosa26^(LSL-Tomato)) (3) a Cre-regulated Cas9 allele (H11^(LSL-Cas9)),as well as (4) homozygous floxed alleles of either p53 or Lkb1.Lentiviral vectors stably integrate into the genome of the transducedcell. Tumors were initiated in KT; Cas9, KPT; Cas9, and KLT; Cas9 miceto generate 31 different genotypes of lung tumors. Mice were analyzedafter 15 weeks of tumor growth. Genomic DNA was extracted from wholelungs, after the addition of barcoded “bench-mark” cell lines, thesgID-barcode region was PCR amplified, deep-sequenced, and analyzed todetermine the relative expansion of each uniquely barcoded tumor usingthe Tuba-seq pipeline.

FIG. 43. Tumor suppression in Kras^(G12D)-driven lung adenocarcinoma invivo. a. Fold change in sgID representation (ΔsgID representation) inKT; Cas9 mice relative to KT mice, which lack Cas9 and therefore shouldnot expand relative to sg Inert. Several sgRNAs (sgIDs) increase inrepresentation, reflecting the increased growth of tumors withinactivation of the targeted tumor suppressor genes. Means and 95%confidence intervals are shown. b,c. The ability to detect tumorsuppressive effects is improved by analyzing individually-barcodedtumors compared to bulk sgRNA representation (ΔsgID representation). (b)Analysis of the relative size of the 95th percentile tumor with eachsgRNA identifies somewhat similar estimates of relative tumor size asbulk ΔsgID representation, which exhibits wider confidence intervals.(c) P-value of the Log-Normal mean (LN mean) measure of relative tumorsize versus P-value ΔsgID representation. Because individual tumor sizesare measured and then properly normalized to eliminate exogenous sourcesof noise, both the 95^(th) percentile and LN Mean metrics identifyfunctional tumor suppressors with greater confidence and precision. p53loss is an exception, as its growth effects are poorly described by aLog-Normal distribution. All P-values are two-sided and obtained via2×10⁶ Bootstrapping permutation tests and a Bonferroni-correction forthe number of investigated tumor suppressors. d-f. Same as in a-c,except for growth effects in KPT; Cas9 mice. Fold change is relative toKT mice, while 95^(th) percentile and LN Mean size estimates arerelative to KPT; Cas9 internal sg Inert controls. g-i. Same as in a-c,except for growth effects in KLT; Cas9 mice. No tumor suppressors wouldhave been identified without Tuba-seq.

FIG. 44. Rb and p53 tumor suppressor cooperativity in lungadenocarcinoma identified by Tuba-seq, confirmed in a mouse model usingCre/lox regulated alleles, and supported by the co-occurrence of RB1 andTP53 mutations in human lung adenocarcinoma. a. Relative LN Mean size ofsg Setd2, sg Lkb1 and sg Rb1 tumors. Rb1 inactivation increase tumorsize less that Setd2 or Lkb1 inactivation in the p53-proficient KT; Cas9background. Conversely, Rb1 inactivation increases tumor size to asimilar extent as Setd2 or Lkb1 inactivation in the p53-deficient KPT;Cas9 background. P-values test null hypothesis of similar LN Mean to sgRb1. P<0.05 in bold. b. H&E staining of representative lung lobes fromKP and KP; Rb1^(flox/flox) mice with tumors initiated withAdeno-CMV/Cre. Mice were analyzed 12 weeks after tumor initiation. Scalebars=500 μm. c. Representative ex vivo μCT images of the lungs from KPand KP; Rb1^(flox/flox) mice are shown. Lung lobes are outlined with adashed white line. d. Quantification of percent tumor area in K;Rb1^(wt/wt), K; Rb1^(flox/flox), KP; Rb1^(wt/wt), and KP;Rb1^(flox/flox) mice. Histological quantification confirms thatRb1-deletion increases tumor burden more dramatically in p53-deficienttumors. * p-value <0.05, n.s.=not significant. Titer of Ad-Cre isindicated. e,f. Co-occurrence of RB1 and TP53 mutations in two humanlung adenocarcinoma genomics datasets: (e) TCGA 2014 dataset, and (f)the GENIE consortium 2017. P-values were calculated using the DISCOVERstatistical independence test for somatic alterations.

FIG. 45. Deep sequencing of targeted genomic loci confirms creation ofindels at all targeted loci and shows selective expansion of cancercells with indels in the strongest tumor suppressor genes. a. Indelabundance in each region targeted by sgRNAs, as determined by deepsequencing of total lung DNA from the targeted regions of four KPT; Cas9mice. Indel abundance is normalized to the median abundance of sgNeo1,sgNeo2, and sgNeo3. Error bars denote range of abundances observed,while dots denote median. Indels were obsevered in all targeted regions.sgp53 is not shown, as its target site is deleted by Cre-mediatedrecombination of the p53^(floxed) alleles. b. Indel abundance asdescribed in (a) versus the 95th percentile tumor size determined byTuba-seq (as described in FIG. 1d ). Each dot represents a single sgRNAin an individual mouse and each mouse is represented by a unique shape.Indel abundance correlated with Tuba-seq size profiles (as expected),however indel abundance does not measure individual tumor sizes andexhibits greater statistical noise. The largest single tumor in thisentire analysis, as determined by Tuba-seq, was an sgCdkn2a tumor thatsimilarly appeared as an outlier in the indel analysis—furthercorroborating faithful analysis of genetic events by Tuba-seq.

FIG. 46. Validation of the redundancy between Setd2 and Lkb1 in mousemodels and in human lung adenocarcinomas. a. Fluorescence dissectingscope images (top) and H&E stained section (bottom) of lung lobes fromKPT and KPT; Cas9 mice with Lenti-sg Setd2#1/Cre or Lenti-sgNeo2/Creinitiated tumors. Mice were analyzed after 9 weeks of tumor growth. Lunglobes are outlined with a white dashed line in fluorescence dissectingscope images. Top scale bars=5 mm. Bottom scale bars=4 mm. b.Quantification of percent tumor area in KPT; Cas9 mice with Lenti-sgSetd2#1/Cre or Lenti-sgNeo2/Cre initiated tumors, and KPT mice withLenti-sg Setd2#1/Cre initiated tumors. Each dot represents a mouse andhorizontal bars are the mean. There is an increase in tumor area betweenKPT; Cas9 and KPT mice with tumors initiated with the same virus, but nodifference between KPT; Cas9 mice tumors initiated with Lenti-sgSetd2#1/Cre and those initiated with Lenti-sgNeo2/Cre, presumably due tohigh mouse-to-mouse variability. Because these lentiviral vectors werebarcoded, we performed Tuba-seq analysis of these mice to quantify thesize of induced tumors. sg Setd2 increased tumor sizes in KPT; Cas9relative to sg Neo2. ** P<0.01, n.s. is not significant. c,d. Same asa,b except for KLT; Cas9 mice with Lenti-sg Setd2#1/Cre orLenti-sgNeo2/Cre initiated tumors. Mice were analyzed after 9 weeks oftumor growth. Top scale bars=5 mm. Bottom scale bars=4 mm. e,f. Theco-occurrence of SETD2 and LKB1 (HGNC name STK11) in two human lungadenocarcinoma genomics datasets: (e) TOGA 2014 dataset¹⁴ (N=229patients), and (f) the GENIE Consortium (N=1563 patients). Two-sidedP-values were calculated using the DISCOVER statistical independencetest.

FIG. 47. Correspondence of Tuba-seq fitness measurements to humangenomic patterns. a. Relative fitness measurements and humanco-occurrence rates of the nineteen pairwise interactions that weinvestigated. LN Mean Ratio is the ratio of relative LN Mean (sg TS/sgInert) within the background of interest divided by the mean relative LNmean of all three backgrounds. Background rate can be either anunweighted average of the three backgrounds (raw), or weighted by eachbackground's rate of occurrence in human lung adenocarcinoma (weighted).*OR=“Odds Ratio” of the co-occurrence rate of a gene pair within thehuman data. One sided P-values of human co-occurrence rates (>0.5suggest mutual exclusivity) were determined using the DISCOVER test.Combined P-values generated using Stouffer's Method (Methods). P<0.025and P >0.975 are in bold. Fitness measurements and co-occurrence ratesgenerally correspond (Spearman's r=0.50, P-value=0.03 for weighted LNMean Ratio; r=0.4 for unweighted). b. Graphical summary of fitnessmeasurements and co-occurrence rates from a. Human GeneticsCooperativity was defined by a Combined Odds Ratio >1 and Redundant <1.Cooperativity for Tuba-seq data denotes a LN Mean Ratio >1 and Redundant<1. c. Number of statistically-significant genetic interactionssuggested from a pan-cancer analysis of twenty-one tumor types. Tumortypes abbreviations are borrowed from TCGA. Lung adenocarcinoma (LUAD)is black and is predicted to contain a quantity of genetic interactionsthat is similar to the median, suggesting that the ruggedness of thefitness landscape studied here may be representative of cancer evolutionin general.

FIG. 48. Power analysis of larger genetic surveys. By assuming lognormaltumor size distributions, the statistical power of Tuba-seq to detectdriver growth effects and non-additive driver interactions in largergenetic surveys can be projected. Future experiments could utilizelarger mouse cohorts and larger pools of sgRNAs targeting putative tumorsuppressors. In all hypothetical experiments, the Lenti-sg TS-Pool/Cretiters and fraction of the pool with inert sgRNAs (for normalization)were kept consistent with our original experiments. a. P-value contoursfor the confidence in detecting a weak driver (parameterized by thesgCdkn2a distribution in KT; Cas9 mice). Any experimental setups above acontour detects weak drivers with a confidence greater than or equal tothe P-value of the contour. b,c. Same as in a, except for moderate andstrong drivers respectively (parameterized by sg Rb1 and sg Lkb1 in KT;Cas9 mice). sgRNA pool size is extended to 500 targets (instead of 100targets in a pool) because larger screens are possible wheninvestigating genes with these effect strengths. d-f. Same as in a-c,except for driver interactions. Driver interactions (LN Mean Ratio) aredefined as a ratio of driver growth rates (sg TS/sg Inert in background#1)/(sg TS/sg Inert in background #2) that were statistically differentfrom the null hypothesis of one. (d) A weak driver interactionparameterized by Rbm10—p53 (7% effect size). (e) A moderate driverinteraction parameterized by Rb1—p53 (13% effect size). (f) A strongdriver interaction parameterized by Setd2—Lkb1 (68% effect size).

DETAILED DESCRIPTION

Before the present methods and compositions are described, it is to beunderstood that this invention is not limited to a particular method orcomposition described, and as such may, of course, vary. It is also tobe understood that the terminology used herein is for the purpose ofdescribing particular embodiments only, and is not intended to belimiting, since the scope of the present invention will be limited onlyby the appended claims.

Where a range of values is provided, it is understood that eachintervening value, to the tenth of the unit of the lower limit unlessthe context clearly dictates otherwise, between the upper and lowerlimits of that range is also specifically disclosed. Each smaller rangebetween any stated value or intervening value in a stated range and anyother stated or intervening value in that stated range is encompassedwithin the invention. The upper and lower limits of these smaller rangesmay independently be included or excluded in the range, and each rangewhere either, neither or both limits are included in the smaller rangesis also encompassed within the invention, subject to any specificallyexcluded limit in the stated range. Where the stated range includes oneor both of the limits, ranges excluding either or both of those includedlimits are also included in the invention.

Unless defined otherwise, all technical and scientific terms used hereinhave the same meaning as commonly understood by one of ordinary skill inthe art to which this invention belongs. Although any methods andmaterials similar or equivalent to those described herein can be used inthe practice or testing of the present invention, some potential andpreferred methods and materials are now described. All publicationsmentioned herein are incorporated herein by reference to disclose anddescribe the methods and/or materials in connection with which thepublications are cited. It is understood that the present disclosuresupersedes any disclosure of an incorporated publication to the extentthere is a contradiction.

As will be apparent to those of skill in the art upon reading thisdisclosure, each of the individual embodiments described and illustratedherein has discrete components and features which may be readilyseparated from or combined with the features of any of the other severalembodiments without departing from the scope or spirit of the presentinvention. Any recited method can be carried out in the order of eventsrecited or in any other order which is logically possible.

It must be noted that as used herein and in the appended claims, thesingular forms “a”, “an”, and “the” include plural referents unless thecontext clearly dictates otherwise. Thus, for example, reference to “acell” includes a plurality of such cells (e.g., a population of suchcells) and reference to “the protein” includes reference to one or moreproteins and equivalents thereof, e.g. polypeptides, known to thoseskilled in the art, and so forth.

The publications discussed herein are provided solely for theirdisclosure prior to the filing date of the present application. Nothingherein is to be construed as an admission that the present invention isnot entitled to antedate such publication by virtue of prior invention.Further, the dates of publication provided may be different from theactual publication dates which may need to be independently confirmed.

Methods and Compositions

As summarized above, compositions and methods are provided for measuringpopulation size for a plurality of clonal cell populations in the sameindividual. As an example, in some cases a subject method is a method ofmeasuring tumor size (e.g., the number of neoplastic cells within atumor) for a plurality of clonally independent tumor cell populations(e.g., different tumors) of the same individual. In some cases a subjectmethod includes: (a) contacting a tissue of an individual with aplurality of cell markers that are heritable and distinguishable fromone another, to generate a plurality of distinguishable lineages ofheritably marked cells within the contacted tissue; (b) after sufficienttime has passed for at least a portion of the heritably marked cells toundergo at least one round of division, detecting and measuringquantities of at least two of the plurality of cell markers present inthe contacted tissue, thereby generating a set of measured values; and(c) calculating, using the set of measured values as input, a number ofheritably marked cells present in the contacted tissue for at least twoof said distinguishable lineages of heritably marked cells.

Contacting a Tissue

In some embodiments, a subject method includes a step of contacting atissue (e.g., a tissue of an individual) (e.g., muscle, lung, bronchus,pancreas, breast, liver, bile duct, gallbladder, kidney, spleen, blood,gut, brain, bone, bladder, prostate, ovary, eye, nose, tongue, mouth,pharynx, larynx, thyroid, fat, esophagus, stomach, small intestine,colon, rectum, adrenal gland, soft tissue, smooth muscle, vasculature,cartilage, lymphatics, prostate, heart, skin, retina, and reproductiveand genital systems, e.g., testicle, reproductive tissue, and the like)with a plurality of cell markers that are heritable and distinguishablefrom one another, to generate a plurality of distinguishable lineages ofheritably marked cells within the contacted tissue. In some cases, thetissue is an engineered tissue grown outside of an animal (e.g., anorganoid, cells in culture, etc.). In some cases, the tissue is part ofa living animal, and therefore the tissue can be considered a tissue ofan individual and said contacting can be performed by administering(e.g., via injection) the cell markers to the individual.

Any convenient route of administration can be used (e.g., intratracheal,intranasal, retrograde pancreatic ductal, intramuscular, intravenous,intraperitoneal, intravesicular, intraarticular, topically,subcutaneous, orally, intratumoral, and the like). In some cases,administration is via injection (e.g., injection of a library, such as aviral library, directly into the target tissue). In some cases, thetransfer of markers into cells is via electroporation (e.g.,nucleofection), transfection (e.g., using calcium phosphate, cationicpolymers, cationic lipids etc), hydrodynamic delivery, sonoporation,biolistic particle delivery, or magnetofection. Any convenient deliveryvector can be used (e.g., viral particles, viral-like particles, nakednucleic acids, plasmids, oligonucelotides, exosomes, lipoplexes,gesicles, polymersomes, polyplexes, dendrimers, nanoparticles, biolisticparticles, ribonucleoprotein complexes, dendrimers, cell-penetratingpeptides, etc.).

The tissue can be any tissue type from any desired animal. For example,in some embodiments the contacted tissue is an invertebrate tissue(e.g., an ectdysozoan, lophotrocozoan, porifera, cnidarian, ctenophoran,arthropod, annelid, mollusca, flatworm, rotifera, arthropod, insect, orworm tissue). In some embodiments the contacted tissue is a vertebratetissue (e.g., an avian, fish, amphibian, reptilian, or mammaliantissue). Suitable tissues also include but are not limited to tissuefrom: rodents (e.g., rat tissue, mouse tissue), ungulates, farm animals,pigs, horses, cows, sheep, non-human primates, and humans. The targettissue can include, but is not limited to: muscle, lung, bronchus,pancreas, breast, liver, bile duct, gallbladder, kidney, spleen, blood,gut, brain, bone, bladder, prostate, ovary, eye, nose, tongue, mouth,pharynx, larynx, thyroid, fat, esophagus, stomach, small intestine,colon, rectum, adrenal gland, soft tissue, smooth muscle, vasculature,cartilage, lymphatics, prostate, heart, skin, retina, and reproductiveand genital systems, e.g., testicle, reproductive tissue, and the like.

In some cases the tissue is contacted for the purpose of inducing cellsto become neoplastic, e.g., in some cases the tissue is contacted forthe purpose of initiating multiple independent tumors to form. Forexample, in some cases the introduced cell markers (and/or componentslinked with the cell markers) cause neoplastic transformation (lead toneoplastic cell formation) and the outcome of multiple differentneoplastic initiating events can be compared to one another because eachevent was uniquely marked with an identifiable heritable cell marker. Insome such cases, the cell markers initiate the same genetic change suchthat the induced tumors begin due to the same type (or even identical)genetic perturbation, but the outcome of each initiating event can betracked because each individual cell marker is distinguishable from theothers. The purpose of such a method may be, for example, to trackmultiple independent cell lineages in the same tissue (and/or sameanimal) in order to generate a population size (e.g., tumor size, numberof neoplastic cells in each tumor) distribution profile for a givengenotype of interest. Alternatively, in some cases different geneticperturbations are used (e.g., the cell makers can cause two or moredifferent genetic perturbations, components linked to the cell makerscan cause two or more different genetic perturbations) and the outcomesfrom different genotypes in the same tissue (e.g., in some cases in thesame animal) can be compared (e.g., different tumors with differentgenetic underpinnings that are present in the same tissue, e.g.,multiple different tumors in the lung, muscle, kidney, and the like).

In some embodiments the tissue already contains neoplastic cells (e.g.,tumors) prior to the contact with the cell markers. In some cases, atumor is contacted with the cell markers (e.g., the cell markers can beinjected into the tumor, injected into the bloodstream to contact thetumor[s], administered to another organ or tissue to contact thetumor[s], etc.). As an example, in some cases, the cell markers are usedas a way to mark independent neoplastic cells such as different cellswithin a neoplasm or tumor, and each marked cell can then be treated asa separate lineage—one can track the number of cells produced for eachtracked lineage by counting the number of cells with each marker present(cells with each marker present) after one or more rounds of celldivision. In some case, the method includes genetically modifying thecells into which the cell markers are introduced. For example, a tissuemay already have one or more tumors prior to performing a subjectmethod, and the purpose of introducing the cell markers is to test theeffect of introducing additional genetic modifications to the tumorcells (i.e., changes in addition to those already present in theneoplastic cells). As such, each distinguishable cell marker can beassociated with a different genetic change (e.g., by pairing nucleicacids encoding guide RNAs that target particular genetic targets with aunique identifier such as a DNA barcodes so that each guide RNA, andtherefore each genetic modification, is associated with a uniqueidentifier such as a DNA barcode). In such a case, the marked lineagerepresents sets of cells that are genetically different (e.g., has amutation at a particular genetic locus) from one another.

Alternatively, in some cases each of the tumors is genetically the sameand the cell markers track lineages that are not necessarily geneticallydifferent from one another. This allows the performer of the method totrack multiple independent cell lineages in the same animal and togenerate a population size (e.g. tumor size, number of neoplastic cellsin tumors) distribution profile for a given genotype of interest.

Cell Markers

A plurality of cell markers (i.e., introduced (heterologous, artificial)cell markers—where the markers are not those that pre-exist in thecells—e.g., the introduced markers are not simply pre-existing clonalsomatic mutations in a tumor) is two or more (e.g., 3 or more, 5 ormore, 10 or more, or 15 or more, 50 or more, 100 or more, 200 or more,500 or more, 1000 or more, 10,000 or more, 100,000 or more, 1,000,000 ormore, 1,000,000,000 or more, etc.) cell markers. Likewise, a pluralityof marked cell lineages is two or more (e.g., 3 or more, 5 or more, 10or more, or 15 or more, 100 or more, 1,000 or more, 10,000 or more,100,000 or more, etc.) marked cell lineages. Any convenient heritablecell markers (that are distinguishable from one another) can be used anda number of heritable cell markers will be known to one of ordinaryskill in the art. In some cases, the cell markers (i.e., introduced(heterologous, artificial) that are heritable and distinguishable fromone another) are barcoded nucleic acids. In some cases, the barcodednucleic acids can be integrated into the genomes of the target cells orin some cases the barcoded nucleic acids can be maintained episomally.Barcoded nucleic acids include nucleotide sequences that provide aunique identifier for each cell lineage that will be detected andquantified/measured. In some case, the plurality of cell markers thatare heritable and distinguishable from one another is a library ofbarcoded nucleic acids, where the exact sequence of the barcode has somerandom element. For example, in some cases the barcode can be describedwith a series of Ns (e.g., positions in the nucleic acid sequence forwhich each nucleotide is not defined but is one of all possible or adefined subset of canonical or non-canonical nucleotides). A subjectbarcoded nucleic acid can include any convenient number of Ns.

In some cases, a subject barcoded nucleic acid (a plurality/library)includes 5 or more (e.g., 6 or more, 7 or more, 8 or more, 10 or more,12 or more, or 15 or more) randomized positions, e.g., 5 or more (e.g.,6 or more, 7 or more, 8 or more, 10 or more, 12 or more, or 15 or more)positions at which the nucleotide is not predetermined. In some casesthe formula for a library (plurality) of barcoded nucleic acids includesa stretch of nucleotides at least 10 base pairs (bp) long (e.g., atleast 12 bp, 15 bp, 17 bp, or 20 bp long) in which 5 or more positions(e.g., 6 or more, 7 or more, 8 or more, 10 or more, 12 or more, or 15 ormore positions) are not defined (i.e., positions at which the baseidentity differs among members of the library). In some cases theformula for a library (plurality) of barcoded nucleic acids includes astretch of nucleotides in which from 5 to 40 positions (e.g., 5 to 30, 5to 25, 5 to 20, 5 to 18, 5 to 15, 5 to 10, 8 to 40, 8 to 30, 8 to 25, 8to 20, 8 to 18, 8 to 15, 8 to 10, 10 to 40, 10 to 30, 10 to 25, 10 to20, 10 to 18, 10 to 15, 12 to 40, 12 to 30, 12 to 25, 12 to 20, 12 to18, or 12 to 15 positions) are not defined (i.e., positions at which thebase identity differs among members of the library). In some cases theformula for a library (plurality) of barcoded nucleic acids includes astretch of nucleotides in which from 5 to 1000 positions (e.g., 5 to800, 5 to 600, 5 to 500, 5 to 250, 5 to 150, 5 to 100, 5 to 50, 5 to 30,5 to 25, 5 to 20, 5 to 18, 5 to 15, 5 to 10, 8 to 1000, 8 to 800, 8 to600, 8 to 500, 8 to 250, 8 to 150, 8 to 100, 8 to 50, 8 to 40, 8 to 30,8 to 25, 8 to 20, 8 to 18, 8 to 15, 8 to 10, 10 to 1000, 10 to 800, 10to 600, 10 to 500, 10 to 250, 10 to 150, 10 to 100, 10 to 50, 10 to 40,10 to 30, 10 to 25, 10 to 20, 10 to 18, 10 to 15, 12 to 1000, 12 to 800,12 to 600, 12 to 500, 12 to 250, 12 to 150, 12 to 100, 12 to 50, 12 to40, 12 to 30, 12 to 25, 12 to 20, 12 to 18, or 12 to 15 positions) arenot defined (i.e., positions at which the base identity differs amongmembers of the library).

The barcoded nucleic acids can be linear (e.g., viral) or circular(e.g., plasmid) DNA molecules. The barcoded nucleic acids can besingle-stranded or double-stranded DNA molecules. Non-limiting examplesinclude plasmids, synthesized nucleic acid fragments, synthesizedoligonucleotides, minicircles, and viral DNA. Barcoded nucleic acids canbe RNA molecules, DNA (DNA molecules), RNA/DNA hybrids, or nucleicacid/protein complexes.

In some cases, cell markers may include a plurality of biomarkers (e.g.,antibodies, fluorescent proteins, cell surface proteins) that areheritable and distinguishable from each other, alone or in combinationwith a plurality of other biomarkers of the same or different type, thatare distinguishable from each other as well as distinguishable from theplurality of other biomarkers when used in combination. In such as case,the biomarkers may be present in a predefined or randomized manner,inside or outside individual cells and/or cell lineages, and can bequantified and/or measured using methods that will be commonly known byone of ordinary skill in the art (e.g. high-throughput/next-generationDNA sequencing, microscopy, flow-cytometry, mass spectrometers, etc).

Cell markers can be delivered to cells using any convenient method. Insome cases, the cell markers (e.g., barcoded nucleic acids) aredelivered to the tissue via viral vector. Any convenient viral vectorcan be used and examples include but are not limited to: lentiviralvectors, adenoviral vectors, adeno-associated viral (AAV) vectors,bocavirus vectors, foamy virus vectors, and retroviral vectors.

In one example from the working examples below (see FIG. 4a ), theplurality of cell markers was delivered to the target tissue vialentiviral vectors. A library of lentiviral particles was used in whicheach viral particle included one barcoded nucleic acid that included atwo-component barcode, where the first component was unique to eachencoded guide RNA and the second component was unique to each moleculeso that in turn it would be unique to each cell lineage that was to bedetected and quantified/measured. The formula for the sequence of thebarcode's second component was NNNNNTTNNNNNAANNNNN. Thus, in a stretchof 19 base pairs, 15 of them were not defined (e.g., randomized). Eachbarcoded nucleic acid of the library: (i) encoded a CRISPR/Cas guideRNA; (ii) included a first barcode—a unique identifier 8-nucleotidebarcode that was linked to the guide RNA such that each different guideRNA sequence was linked to its own unique 8-nucleotide barcode; (iii)included a second barcode—the random 19 nucleotide barcode above with 15undefined positions [for tracking cell lineage]; and (iv) encoded a geneediting protein (CRE), the expression of which would lead to Cas9expression in the target tissue. Thus, in this case, multiple differentmembers of the plurality of barcoded nucleic acids included the samefirst barcode, where each first barcode had a ‘corresponding’ guide RNA.However, the second barcode was unique to each member of the librarysuch that each cell lineage that will be detected andquantified/measured would have a unique identifier. Thus, while somemembers shared a first barcode sequence because they shared a commonguide RNA, each member of the library had a unique second barcode thatcould be used to track each integration (i.e., each lineage).

In some cases, a plurality of cell markers that are heritable anddistinguishable are associated with one or more (e.g., 1 or more, 2 ormore, 3 or more, 5 or more, 7 or more, 9 or more, 11 or more, 13 ormore, 15 or more, or 20 or more) pluralities of cell markers that areheritable and distinguishable from one another as well asdistinguishable from the cell markers of the other pluralities of cellmarkers they are associated with. For example, one barcoded nucleic acidmay include a four-component barcode, where the first component isunique to a candidate therapy (e.g. candidate anti-cancer compound), thesecond component is unique to each individual (e.g. a mouse who may ormay not receive the candidate therapy), the third component is unique toan encoded guide RNA, and the fourth component is unique to eachmolecule, so that in turn, the barcoded nucleic acid would be unique toeach cell lineage that was to be detected and quantified/measured. Thus,in this example, the number of cells in each cell lineage can bequantified/measured and each cell lineage can also be directly linked byits four-component nucleic acid barcode to the specific geneticperturbation induced by the guide RNA in that cell lineage, the specificcandidate therapy encountered by that cell lineage, and the specificindividual (e.g., mouse) within which the cell lineage resided.

In some cases, the barcode is incorporated into a DNA donor template forhomology directed repair (HDR) or, e.g., any other mechanism thatincorporates a defined nucleic acid sequence into a desired position inthe genome. For example, the HDR repair template may be used tointroduce the same coding change (e.g. same coding allele), or even asubset of desired changes, into the genome of the cells it contacts, buteach integration event can be independently tagged because the libraryof HDR templates has been randomized at particular positions. In oneexample from the working examples below (see, e.g., FIG. 23a ), theplurality of cell markers (a library of AAV particles in which each AAVparticle included one HDR template) was delivered to the target tissueby AAV particles. The HDR template in each AAV included one of the 12possible non-synonymous, single-nucleotide point mutations in Krascodons 12 and 13 or the wild type Kras sequence as well as a random8-nucleotide barcode in the wobble positions of the adjacent codons touniquely tag each cell that undergoes HDR. The barcode was(N)GG(N)AA(R)TC(N)GC(N)CT(N)AC(N)AT(H) (SEQ ID NO: 1), and thus was astretch of 22 base pairs in which 8 positions were not defined.

In some cases, the cell markers may contact the tissue in response toexternal perturbation (e.g., candidate anti-cancer therapy). In such acase, the administration of the external perturbagen may occurstochastically, with tunable probabilities, or as a result of acombinatorial matching of signals (e.g., a predefined physiologicalstate of the cell, the level of expression of a specific gene, set ofgenes, or sets of genes, the level of activity of a specific pathway orpathways, and/or other signals internal or external to the cell or celllineage [e.g., the identity of the tissue, levels of blood supply,immune state of the whole individual, physical location of the cell,etc]). For example, a cell marker (e.g. barcoded DNA) may contact thetissue upon expression of both a guide RNA, under the control of anenhancer specific to particular type of epithelial cell, and Cas9, inresponse to a compound being administered to the individual within whichthe tissue exists.

In some cases, the cell markers may contact a healthy or diseased cellpopulation or tissue in vivo in an individual living organism, or invitro in a cell population in culture or an organoid culture. In somecases, cell markers may contact a neoplastic cell lineage that isincreasing or decreasing in number or static. In some cases, cellmarkers may contact the tissue in response to administration of a drugor other physiological or environmental perturbation, stochasticallywith tunable probabilities, or via a counting mechanism that induced thecell marker to contact the tissue after a certain number of celldivisions, exactly or stochastically, with tunable mean and variance andother moments, or as a result of a combinatorial matching of signals.

Genetic Modification (Alteration) of Target Cells

As noted above, in some embodiments, the method includes geneticallymodifying the cells into which the cell markers are introduced. In somesuch cases, the introduced cell markers are the agents of the geneticmodification. For example, in some cases the cell markers are barcodednucleic acids that induce genetic modification (e.g., genomicmodification) and in some such cases are barcoded nucleic acids thatinduce neoplastic cell formation. For example, expression of an RNA(e.g., guide RNA) and/or protein (e.g., Cre, a CRISPR/Cas RNA-guidedprotein, etc.) from the barcoded nucleic acids can lead to one or moregenomic alterations, and in some cases the genomic alterations result intransformation of the target cell into a neoplastic cell (e.g., which insome cases can result in tumor formation).

However, whether a cell marker (e.g., barcoded nucleic acid) introducesa genomic modification can be independent of whether it can induceneoplastic cell formation. For example, in some cases a barcoded nucleicacid can encode an oncogene (a gene that when expressed as a protein canlead to neoplastic cell formation). In some such cases, the barcodednucleic acid does not induce a genomic change in the target cell butdoes induce neoplastic cell formation due to expression of the oncogene.In some cases, an oncogene encodes a wild type protein that can cause acell to become neoplastic when the protein is overexpressed. In somecases an oncogene encodes a mutated protein (e.g., mutated form of KRAS)that can cause a cell to become neoplastic when the protein isexpressed. In some cases a cell marker (e.g., barcoded nucleic acid)introduces a genomic modification in the target cell but themodification only induces neoplastic formation (e.g., tumor/cancerformation) in combination with one or more additional genomicmodifications that may occur before, during, or a period time after theintroduction of the cell marker and associated genomic modification.

On the other hand, in some cases a cell marker (e.g., barcoded nucleicacid) introduces a genomic modification in the target cell but themodification does not induce neoplastic formation (e.g., tumor/cancerformation). For example in some cases a barcoded nucleic acid integratesinto the genome of a target cell in an inert way.

In some cases a barcoded nucleic acid encodes a protein (e.g., wild typeor mutant protein) where the protein is not necessarily related tocancer, e.g., the protein(s) can be involved in any biological processof interest and its expression may not have an effect on cellproliferation and/or neoplastic cell formation (e.g., may not be anoncogene or a tumor suppressor). In some such cases the nucleic acidintegrates into the genome of target cells and in other cases thenucleic acid does not integrate into the genome (e.g., can be maintainedepisomally). In some cases a barcoded nucleic acid encodes wild type ormutant protein, e.g., a cDNA, that encodes a protein that is detrimentalto tumors, e.g., in some way other than growth/proliferation control.

In some embodiments a subject cell marker (e.g., barcoded nucleic acid)both introduces a genomic modification in the target cell and alsoinduces neoplastic cell formation (e.g., tumor/cancer formation). Forexample, in some cases a barcoded nucleic acid can cause editing at atarget locus to modify a tumor suppressor, alter the expression of anoncogene, edit a gene (e.g., Kras) to become a neoplastic-inducingallele, etc.

As noted above, expression of an RNA (e.g., guide RNA) and/or protein(e.g., Cre, a CRISPR/Cas RNA-guided protein, etc.) from a barcodednucleic acid can lead to one or more genomic alterations, and in somecases the genomic alterations result in transformation of the targetcell into a neoplastic cell (e.g., which in some cases can result intumor formation). In some embodiments, genomic alteration of the targetcells can be temporally separated from the initiation of neoplasticcharacter (e.g., from tumor initiation). As an example, a vector(s)could be engineered to allow temporal control of a CRISPR/Cas guide RNAand/or temporal control of CRISPR/Cas nucleic acid-guided proteinactivity (e.g., Cas9 activity).

In some cases, a protein that introduces genetic (e.g., genomic)modification is expressed in the target cells. The protein can beintroduced into a target cell as protein or as a nucleic acid (RNA orDNA) encoding the protein. The protein may also already be encoded by anucleic acid in the cell (e.g., encoded by genomic DNA in the cell) andthe method includes inducing the expression of the protein. In somecases a protein that introduces a genetic modification in target cellsof a target tissue is a genome editing protein/endonuclease (some ofwhich are ‘programmable’ and some of which are not). Examples includebut are not limited to: programmable gene editing proteins (e.g.,transcription activator-like (TAL) effectors (TALEs), TALE nucleases(TALENs), zinc-finger proteins (ZFPs), zinc-finger nucleases (ZFNs),DNA-guided polypeptides such as Natronobacterium gregoryi Argonaute(NgAgo), CRISPR/Cas RNA-guided proteins such as Cas9, CasX, CasY, Cpf1,and the like) (see, e.g., Shmakov et al., Nat Rev Microbiol. 2017 March;15(3):169-182; and Burstein et al., Nature. 2017 Feb. 9;542(7640):237-241); transposons (e.g., a Class I or Class IItransposon—e.g., piggybac, sleeping beauty, Tc1/mariner, Tol2,PIF/harbinger, hAT, mutator, merlin, transib, helitron, maverick, frogprince, minos, Himarl and the like); meganucleases (e.g., I-SceI,I-CeuI, I-CreI, I-DmoI, I-ChuI, I-Dirt, I-FlmuI, I-FlmuII, I-AniI,I-SceIV, I-CsmI, I-Pant, I-PanII, I-PanMI, I-SceII, I-PpoI, I-SceIII,I-LtrI, I-GpiI, I-GZeI, I-OnuI, I-HjeMI, I-MsoI, I-TevI, I-TevII,I-TeviII, PI-MleI, PI-MtuI, PI-PspI, PI-Tli I, PI-Tli II, PI-SceV, andthe like); megaTALs (see, e.g., Boissel et al., Nucleic Acids Res. 2014February; 42(4): 2591-2601); phage-derived integrases; a Cre protein; aFlp protein; and the like. In some cases the genome editing nuclease(e.g., a CRISPR/Cas RNA-guided protein) has one or more mutations thatremove nuclease activity (is a nuclease dead protein) and the protein isfused to a transcriptional activator or repressor polypeptide (e.g.,CRISPRa/CRISPRO. In some cases the genome editing nuclease (e.g., aCRISPR/Cas RNA-guided protein) has one or more mutations that removenuclease activity (is a nuclease dead protein) or partially removenuclease activity (is a nickase protein), may have one or moreadditional mutations that modulate protein function or activity, and theprotein is fused to a deaminase domain (e.g., ADAR, APOBEC1, etc.),which itself may have one or more additional mutations that modulateprotein function or activity, or fused to the deaminase domain and oneor more additional proteins or peptides (e.g., the bacteriophage Gamprotein, uracil glycosylase inhibitor, etc.), which may also have one ormore additional mutations that modulate protein function or activity(e.g., RNA base editors, DNA base editors).

In some cases, an editing protein such as Cre or Flp can be introducedinto the target tissue for the purpose of inducing expression of anotherprotein (e.g., a CRISPR/Cas RNA-guided protein such as Cas9) from thegenome, e.g., an animals can contain a lox-stop-lox allele of Cas9 andan introduced Cre protein (e.g., encoded by a barcoded nucleic acid)results in removal of the ‘stop’ and thus results in expression of theCas9 protein.

In some embodiments, the barcoded nucleic acids can induce neoplasticcell formation and include one or more of: homology directed repair(HDR) DNA donor templates, nucleic acids encoding oncogenes (includingwild type and/or mutant alleles of proteins), nucleic acids encodingCRISPR/Cas guide RNAs, nucleic acids encoding short hairpin RNAs(shRNAs), and nucleic acids encoding a genome editing protein (e.g., seeabove).

In some cases when the barcoded nucleic acids are HDR DNA donortemplates, they can introduce mutations into the genome of target cells.In some such cases, a genome editing nuclease is present in the cell(either introduced or induces as part of the subject method or alreadyexpressed in the targeted cells) that will cleave the targeted DNA suchthat the donor templates are used to insert the barcoded sequence. Insome cases, a library (plurality) of HDR DNA donor templates includesmembers that have unique sequence identifiers (barcodes) for eachmolecule, but the molecules result in the same functional perturbation(e.g., they may all result in expression of the same protein, e.g., insome cases with a mutated amino acid sequence, but they may differ inthe wobble positions of the codons then encode the protein such that theresulting multiple cell lineages are distinguishable from one anotherdespite expressing the same mutated protein). In some cases, a library(plurality) of HDR DNA donor templates includes members that have uniquesequence identifiers (barcodes) for each molecule, and the moleculesresult in the different functional perturbations (e.g., can targetdifferent genetic loci, can target the same loci but introduce differentalleles, etc.).

In some cases the barcoded nucleic acids are CRISPR/Cas guide RNAs orare DNA molecules that encode CRISPR/Cas guide RNAs. A library of suchmolecules can include molecules that target different loci and/ormolecules that target the same locus. In some cases the barcoded nucleicacids encode an oncogene, which for purposes of this disclosure includeswild type proteins that can cause neoplastic cell formation whenoverexpressed as well as mutated proteins (e.g., KRAS—see workingexamples below) that can cause neoplastic cell formation. A library ofsuch molecules can include molecules that express the same oncogene or alibrary of molecules that express different oncogenes. In some cases thebarcoded nucleic acids include short hairpin RNAs (shRNAs) and/or DNAmolecule(s) that encode shRNAs (e.g., which can be targeted to anydesired gene, e.g., tumor suppressors). A library of such molecules caninclude molecules that express the same shRNAs or a library of moleculesthat express different shRNAs. In some cases the barcoded nucleic acidsinclude RNAs and/or DNAs that encode one or more genome editingproteins/endonucleases (see above for examples, e.g., CRISPR/CasRNA-guided proteins such as Cas9, Cpf1, CasX or CasY; Cre recombinase;Flp recombinase; ZFNs; TALENs; and the like). A library of suchmolecules can include molecules that express the same genome editingproteins/endonucleases or a library of molecules that express differentgenome editing proteins/endonucleases.

In some embodiments, the cell markers are distinguishably labeledparticles (e.g., beads, nanoparticles, and the like). For example, insome cases the particles can be labeled with distinguishable mass tags(which can be analyzed via mass spectrometry), with distinguishablefluorescent proteins, with distinguishable radio tags, and the like.

Detecting/Measuring/Calculating

Subject methods can also include, e.g., after sufficient time has passedfor at least a portion of the heritably marked cells to undergo at leastone round of division, a step of detecting and measuring quantities ofat least two of the plurality of cell markers present in the contactedtissue.

In some cases, the period time that elapsed between steps (a) and (b)[between contacting a tissue with a plurality of cell makers anddetecting/measuring the cell markers present in the tissue] is a periodof time sufficient for at least a portion (e.g., at least two of thedistinguishably marked cells) of the heritably marked cells to undergoat least one round of division (e.g., at least 2 rounds, 4 rounds, 6rounds, 8 rounds, 10 rounds, or 15 rounds of cell division). In somecases, the period time that elapsed between steps (a) and (b) [betweencontacting a tissue with a plurality of cell makers anddetecting/measuring the cell markers present in the tissue] is 2 or morehours (e.g., 4 or more, 6 or more, 8 or more, 10 or more, 12 or more, 15or more, 18 or more, 24 or more, or 36 or more hours). In some cases,the period time that elapsed between steps (a) and (b) [betweencontacting a tissue with a plurality of cell makers anddetecting/measuring the cell markers present in the tissue] is 1 or moredays (e.g., 2 or more, 3 or more, 4 or more, 5 or more, 7 or more, 10 ormore, or 15 or more, 20 or more, or 24 or more days). In some cases, theperiod time that elapsed between steps (a) and (b) [between contacting atissue with a plurality of cell makers and detecting/measuring the cellmarkers present in the tissue] is 1 or more week (e.g., 2 or more, 3 ormore, 4 or more, 5 or more, 7 or more, or 10 or more weeks). In somecases, the period time that elapsed between steps (a) and (b) [betweencontacting a tissue with a plurality of cell makers anddetecting/measuring the cell markers present in the tissue] is in arange of from 2 hours to 60 weeks (e.g., from 2 hours to 40 weeks, 2hours to 30 weeks, 2 hours to 20 weeks, 2 hours to 15 weeks, 10 hours to60 weeks, 10 hours to 40 weeks, 10 hours to 30 weeks, 10 hours to 20weeks, 10 hours to 15 weeks, 18 hours to 60 weeks, 18 hours to 40 weeks,18 hours to 30 weeks, 18 hours to 20 weeks, 18 hours to 15 weeks, 1 dayto 60 weeks, 1 day to 40 weeks, 1 day to 30 weeks, 1 day to 20 weeks, 1day to 15 weeks, 3 days to 60 weeks, 3 days to 40 weeks, 3 days to 30weeks, 3 days to 20 weeks, 3 days to 15 weeks, 1 week to 60 weeks, 1week to 40 weeks, 1 week to 30 weeks, 1 week to 20 weeks, or 1 week to15 weeks). In some cases, the period time that elapsed between steps (a)and (b) [between contacting a tissue with a plurality of cell makers anddetecting/measuring the cell markers present in the tissue] is in arange of from 2 hours to 300 weeks (e.g., from 2 hours to 250 weeks, 2hours to 200 weeks, 2 hours to 150 weeks, 2 hours to 100 weeks, 2 hoursto 60 weeks, 2 hours to 40 weeks, 2 hours to 30 weeks, 2 hours to 20weeks, 2 hours to 15 weeks, 10 hours to 300 weeks, 10 hours to 250weeks, 10 hours to 200 weeks, 10 hours to 150 weeks, 10 hours to 100weeks, 10 hours to 60 weeks, 10 hours to 40 weeks, 10 hours to 30 weeks,10 hours to 20 weeks, 10 hours to 15 weeks, 18 hours to 300 weeks, 18hours to 250 weeks, 18 hours to 200 weeks, 18 hours to 150 weeks, 18hours to 100 weeks, 18 hours to 60 weeks, 18 hours to 40 weeks, 18 hoursto 30 weeks, 18 hours to 20 weeks, 18 hours to 15 weeks, 1 day to 300weeks, 1 day to 250 weeks, 1 day to 200 weeks, 1 day to 150 weeks, 1 dayto 100 weeks, 1 day to 60 weeks, 1 day to 40 weeks, 1 day to 30 weeks, 1day to 20 weeks, 1 day to 15 weeks, 3 days to 300 weeks, 3 days to 250weeks, 3 days to 200 weeks, 3 days to 150 weeks, 3 days to 100 weeks, 3days to 60 weeks, 3 days to 40 weeks, 3 days to 30 weeks, 3 days to 20weeks, 3 days to 15 weeks, 1 week to 300 weeks, 1 week to 250 weeks, 1week to 200 weeks, 1 week to 150 weeks, 1 week to 100 weeks, 1 week to60 weeks, 1 week to 40 weeks, 1 week to 30 weeks, 1 week to 20 weeks, or1 week to 15 weeks).

The amount (level) of signal detected for each distinguishable cellmarker (e.g., barcoded nucleic acid) can be used to determine the numberof cells present in the contacted tissue (the tissue into which theheritable cell markers were introduced). Any convenient method can beused to detect/measure the cell markers, and one of ordinary skill inthe art will understand that the type of cell markers used will drivewhat method should be used for measuring. For example, if mass tags areused, then mass spectrometry may be the method of choice for measuring.If barcoded nucleic acids are used as the cell markers, then sequencing(e.g., high-throughput/next generation sequencing) may be the method ofchoice for measuring. In some cases, high-throughput sequencing is usedand the number of sequence reads for each detected barcode can be usedto determine the number of cells that contained that particular barcode.In some case the metric of importance is not the number of cells in eachlineage but rather the number of clonal lineages that exceed a certainnumber of cells.

In some cases, sequencing (e.g., high-throughput/next generationsequencing) is performed on PCR products, where the PCR products arefrom PCR reactions that amplified the barcode region from the cellmarkers within the cells (in some cases from the genomic region in whichbarcoded nucleic acids integrated) (see, e.g., FIG. 1a ).

In some cases, the quantification of the number of neoplastic cells intumors, as well as additional phenotyping and analysis, is conductedfrom pooled samples, samples sorted via single, multiple, orcombinatorially arranged biomarkers (e.g., fluorescent proteins,cell-surface proteins, and antibodies), or via dissection of individualtumors from the tissue, organ, cell culture, or other possible means ofcell propagation.

In some cases, ‘benchmarks’ can be used to aid in calculating a cellnumber. For example, in some cases controls can be ‘spiked’ into thesample. For example, spiked (spike in) controls can be used to determinethe number of sequence reads per cell (e.g., number of cells persequence read). For example, in some cases a spiked (spike in) controlcan also be used to correlate the amount of measured DNA with the numberof cells from which the DNA was derived. For example, a known number ofcells can be used to prepare DNA, and the DNA can be processed inparallel with DNA extracted from cells of the contacted tissue (tissuecontacted with heritable cell markers according the methods of thedisclosure). Such a spiked (spike in) control (a ‘benchmark’) caninclude its own unique barcode. The results from the spiked controls canbe used to derive/calculate the number of cells represented by thenumber of sequence reads detected in the sequencing reaction (i.e.,spiked (spike in) controls can be used to provide a coefficient forconverting amount of measured value, e.g., number of sequence reads,into a cell number, e.g., an absolute cell number). Such a process canbe referred to as ‘normalizing’, e.g., sequencing results provide anumber of reads for each unique barcode that is detected, and this valuecan then be compared to one or more ‘benchmarks’ in order calculate anabsolute number of cells that had included the detected unique barcodes(see, e.g., FIG. 1a ).

In some cases because multiple clonal cell populations are detectable bycontacting a subject tissue with the heritable cell markers, and in somecases each distinguishable cell population has a similar genotype, thesubject methods can be used to provide a distribution of population size(e.g., a distribution of tumor size) for a particular phenotype. Forexample, if the initial contacting causes a similar genomic alterationin all contacted cells (e.g., if all cells receive a guide RNA targetingthe same locus, if all cells receive a nucleic acid encoding the sameoncogene allele, and the like), but each cell population (e.g., tumor)is independent, the resulting cell population sizes can provide a clonalcell population size distribution for that particular genotype. Forexample, the goal of performing a subject method may be to search forgenetic changes that alter tumor behavior in particular ways (e.g.,change the size distribution without change the number of tumors perse). For example, the working examples below (e.g., see workingexample 1) include a demonstration that animals with tumors havingp53-deficiency generated a tumor size distribution that was power-lawdistributed for the largest tumors (consistent with a Markov processwhere very large tumors are generated by additional, rarely acquireddriver mutations). Conversely, animals with tumors having Lkb1inactivation increased the size of a majority of lesions suggesting anordinary exponential growth process (e.g., see FIGS. 10, 13, 16, and20).

Size distribution measurements can be used in a number of differentways. For example, one can determine the baseline size distribution ofcell population size (e.g., tumor size) for a given genotype byperforming the methods described herein, and compare it to the sizedistribution that is measured when similarly treated animals are alsotreated with a test compound (e.g., candidate anti-cancer therapy). Thechange in size distribution can be used as a measure of whether the testcompound was effective. As an illustrative example, the inventorsdetermined a baseline measurement for tumor size distribution for micewith tumors that had p53-deficiency, and found that p53-deficiencytended to lead to some tumors that were much larger compared to othertumors. Thus, the size distribution of the p53-deficient tumors was nota standard distribution but instead included outlier tumors. Using themethods described herein, it is possible to screen for potentialtherapeutics (e.g., small molecules, large molecules, radiation, chemo,fasting, antibodies, immune cell therapies, enzymes, viruses, biologics,compounds, and the like) that change the tumor size distribution, but donot necessarily cure animals of tumors. For example, a therapy (e.g., acompound) can be found that, although it does not eradicate allp53-deficient tumors, it instead inhibits the outlier large tumors fromforming. Such a change may not be detected using standard methodsbecause the tested compound would not necessarily reduce overall tumornumber (tumor burden) or even average tumor size (and such a compoundmight be discarded using other methods as a compound that has no effecton inhibiting tumor growth)—but such a therapy (e.g., compound) may bevery useful in clinical settings to treat patients with p53-deficienttumors because it would be effective against the most advanced tumors(e.g., the biggest, more dangerous tumors)(e.g., reduce the risk ofoutlier tumors).

Thus, in some cases, subject methods can be used for screening candidatetherapies (e.g., small molecules, large molecules, radiation,chemotherapy, fasting, antibodies, immune cell therapies, enzymes,viruses, biologics, compounds, and the like) for their effect onpopulation size (e.g., the growth/proliferation of tumors). For example,a subject method can be performed in the presence of a test therapy,e.g., compound (e.g., drug)(e.g., the method can include a step ofcontacting the tissue, e.g., via administration to an individual, withthe test compound), and the effect of the drug can be measured, e.g.,via comparison to parallel experiments in which no drug (e.g., controlvehicle) was added. In cases where the lineage marked cell populationsare genetically the same (or similar) such a method can test whether thecompound has an effect on size distribution of the cell populations. Incases where the differentially marked cells have different genotypes(e.g., different genes have been mutated and/or are being expressed inthe different cell lineages), the therapy (e.g., compound)can be testedagainst multiple different genotypes at the same time, e.g., in the sameanimal in cases where the tissue is in a living animal in vivo. In somecases, such experiments and/or therapy (e.g., compound)screens can beperformed on tissues grown in culture (e.g., 2D cultured tissue, 3Dcultured tissue, organoid cultures). In some cases, such methods can beperformed in non-human animals such as rodents (e.g., mice, rats), pigs,guinea pigs, non-human primates, and the like.

Any perturbagen (e.g., small molecules, large molecules [e.g. antibodiesor decoy receptors], radiotherapies, chemotherapies, inducers ofinflammation, hormones, nanoparticles, immune cell therapies, enzymes,viruses, environmental interventions (e.g. intermittent fasting, acuteexercise, diet control), and the like) can be assessed for its effect onpopulation size for a plurality of marked cell populations. Geneticperturbations can also be induced in all clonal lineages to assess theirimpact. In the case where all lineages are of the same initial genotype,then the response of individual clonal lineages (e.g. tumors) can bedetermined. In the case where the clonal lineages have been induced tohave different defined alterations, then the impact of the induciblegenetic perturbations on clonal lineages with different alterations canbe determined. Systems to generate inducible genetic alteration includebut are not limited to the use of the Flp/FRT or Cre/loxP systems (incell lineages that have not been initiated with Flp or Cre-regulatedalleles) or tetracycline regulatable systems (e.g. tTA or rtTA withTRE-cDNA(s) and/or TRE-shRNA(s) and/or TRE-sgRNA(s)). RegulatableCRISPR/Cas9 genome editing and secondary transduction of neoplasticcells could generated genomic alterations in a temporal manner.

In some cases, the effect of and response to (e.g. pharmacological,chemical, metabolic, pharmacokinetic, immunogenic, toxicologic,behavioral, etc.) an external perturbagen (e.g. candidate anti-cancertherapy) by an individual with a plurality of marked cell populationswill be assessed before, during, and/or after the measuring of cellmarkers.

In some embodiments, a subject method includes, after generatingheritably marked cells (e.g., heritably marked tumors), transplantingone or more of the marked cell populations (e.g., all or part of a tumoror tumors) into a recipient (e.g., a secondary recipient) or a pluralityof recipients, e.g., to seed tumors in the recipient(s). In some cases,such a step can be considered akin to ‘replica plating,’ where one canscreen a large number animals against a test compound, where each animalis seeded from cells from the same starting tumor. Thus, in some cases,the method includes a step in which a test compound is administered tothe recipient(s) of the transplant (e.g., the method can includedetecting and measuring quantities of at least two of the plurality ofcell markers present in the secondary recipient), e.g., to assess growthof the transplanted cells (and some cases this can be done in thepresence and/or absence of a test compound). Thus, a subject method canbe used as part of serial transplantation studies, where the initiallygenerated heritably marked cells (e.g., heritably marked tumors) aretransplanted into one or more recipients, and the number of heritablymarked cells present in the contacted tissue can be calculated for atleast two of the distinguishable lineages of heritably marked cells. Insome of the above cases (e.g., serial transplant) a test compound can beadministered to the serial transplant recipient and the results can becompared to controls (e.g., animals that received a transplant but notthe test compound, animals that received test compound but nottransplant, and the like).

In some embodiments, one or more heritably marked cells are re-marked(e.g., re-barcoded). In other words, in some cases, a population ofcells (e.g., a tumor) that has already been heritably marked iscontacted with a second plurality of cell markers that are heritable anddistinguishable from one another as well as distinguishable from thecell markers of the first plurality of cell markers. In this way, theuser can investigate, e.g., the variability present within marked cellpopulations (e.g., tumors). In some embodiment the heritable markeritself changes over time to record the phylogeny of the cells with aclonal lineage (e.g. evolving nucleotide barcodes).

The heritable lineage marker can also be encoded within an expressedgene (either endogenous or engineered) which facilitates the celllineage to be determined through analysis of mRNA or cDNA from themarked cells. In some cases, cell markers are converted into a differenttype of cell markers (e.g. barcoded DNA expressed by a marked cell asbarcoded RNA or protein). In such a case, one of ordinary skill in theart will understand that the method used for measuring the cell markerwill be determined by the type of cell marker desired to be measured, atthe time of measurement. For example, if barcoded DNA is used as thecell marker and the barcoded DNA is expressed as barcoded RNA, then RNAsequencing (e.g., whole transcriptome sequencing, single cell RNAsequencing, etc.) may be the method used for measuring if RNA barcodesare the type of cell marker that are desired to be measured, or DNAsequencing (e.g., whole genome sequencing, whole exome sequencing,targeted DNA sequencing, etc.) may be the method used for measuring ifDNA barcodes are the type of cell marker that are desired to bemeasured. In such a case, the choice of cell marker to measure may bedriven by the desired phenotype of the cells to investigate and directlylink to cell markers (e.g. barcoded RNA cell markers may be measuredusing single cell RNA sequencing so the RNA expression pattern can bedirectly linked to the cell marker). In some cases, cell lineage markerscan be measured using single cell analysis methods (e.g. single cellRNA-seq, flow cytometry, mass cytometry (CyTOF), MERFISH, single cellproteomics) such that individual cells from each lineage can be relatedto individual cells from each other lineage. In such a case, thephenotypes of the cells within each lineage are investigated. In such acase, these analyses can also used to assess the phenotypic response ofcells of different lineages to external perturbations (e.g., drugtreatment).

When detecting and measuring a heritable cell marker (e.g., a barcodednucleic acid), in some cases the measurement is derived from a wholetissue. As such, a tissue sample can be a portion taken from a tissue,or can be the entire tissue (e.g., a whole lung, kidney, spleen, blood,pancreas, etc.). As such, cell markers (e.g., nucleic acids) can beextracted from a tissue sample so as to represent the remaining tissueor can be extracted from an entire tissue.

In some cases, a biological sample is a blood sample. In some cases thebiological sample is a blood sample but the contacted tissue was not theblood. For example, in some cases a heritably marked cell can secrete acompound (e.g. a unique secreted marker such as a protein or nucleicacid) into the blood and the amount of the compound present in the bloodcan be used to calculate the number of cells present that secret thatparticular compound. For example, heritably marked cells can in somecases secret a fluorescent protein into the blood, and the fluorescentprotein can be detected and measured, and used to calculate the cellpopulation size for cells secreting that particular compound. In somecases, these secreted heritable markers are detected in unperturbedindividuals or after administration of an external perturbagen (e.g.drug).

In some cases, a biological sample is a bodily fluid (e.g., blood, bloodplasma, blood serum, urine, saliva, fluid from the peritoneal cavity,fluid from the pleural cavity, cerebrospinal fluid, etc.). In some casesthe biological sample is a bodily fluid but the contacted tissue was notthe bodily fluid. For example, in some cases a heritably marked cell canrelease an analyte (e.g. a unique marker such as a protein, nucleicacid, or metabolite) into the urine and the amount of the compoundpresent in the urine can be used to calculate the number of cells ornumber of cell lineages that released that particular compound, eitherin alone or in response to an external perturbagen (e.g. candidateanti-cancer therapy).

In some cases, the measuring of cell markers in a biological sample isperformed in parallel with the analysis of cells, cellular components(e.g. cell-free DNA, RNA, proteins, metabolites, etc.), or any otheranalytes (e.g. DNA, RNA, proteins, metabolites, hormones, dissolvedoxygen, dissolved carbon dioxide, vitamin D, glucose, insulin,temperature, pH, sodium, potassium, chloride, calcium, cholesterol, redblood cells, hematocrit, hemoglobin, etc.) that may be directly orindirectly associated with the cell markers and that may be present inthe same biological sample or in a separate biological sample.

In some cases, as noted above, the detecting and measuring is performedon a biological sample collected from an individual (e.g., a bloodsample). In some cases, the detecting and measuring is performed on atissue sample of the contacted tissue, which can in some cases be aportion of the contacted tissue or can be the whole tissue.

When detecting and measuring, biomarkers (other than the introducedheritable cell markers) can be taken in to account. For example, asubject method can include a step of detecting and/or measuring abiomarker of the heritably marked cells, and categorizing the heritablymarked cells based on the results of the biomarker measurements. Such abiomarker can indicate any of number of cellular features, e.g.,proliferation status (e.g., detection of Ki-67 protein, BrdUincorporation, etc.), cell type (e.g., using biomarkers of various celltypes), developmental cell lineage, stemness (e.g., whether a cell is astem cell and/or what type of stem cell), cell death (e.g. Annexin Vstaining, cleaved caspase 3, TUNEL, etc), and cellular signaling state(e.g., detecting phosphorylation state of signaling proteins, e.g.,using phospho-specific antibodies).

In some cases understanding the genotype specificity of a certaintherapy or perturbation can be used to inform (by similarity to othertherapies or perturbations) the mechanism of action of that therapy orperturbation. By uncovering the genotype specificity the methodsdisclosed herein can be used to make and test prediction of combinationtherapies for defined genotypes. Panels of therapies can be tested toestablish their genotype specifity.

Kits and Systems

Also provided are kits and systems, e.g., for practicing any of theabove methods. The contents of the subject kits and/or systems may varygreatly. A kit and/or system can include, for example, one or more of:(i) a library of heritable cell markers that are distinguishable fromone another (e.g., barcoded nucleic acids); (ii) directions forperforming a subject method; (iii) software for calculating the numberof cells from values generated from the detecting and measuring steps ofthe subject methods; (iv) a computer system configured.

In addition to the above components, the subject kits can furtherinclude instructions for practicing the subject methods. Theseinstructions may be present in the subject kits in a variety of forms,one or more of which may be present in the kit. One form in which theseinstructions may be present is as printed information on a suitablemedium or substrate, e.g., a piece or pieces of paper on which theinformation is printed, in the packaging of the kit, in a packageinsert, etc. Yet another means would be a computer readable medium,e.g., diskette, CD, flash drive, etc., on which the information has beenrecorded. Yet another means that may be present is a website addresswhich may be used via the internet to access the information at aremoved site. Any convenient means may be present in the kits.

Utilities

Examples of various applications of the subject matter of thisdisclosure include, but are not limited to the following:

Quantifying the Effect of More Complex Genotypes:

The inventors have already generated and validated lentiviral vectorsthat express pairs of CRISPR/Cas single guide RNAs (sgRNAs),facilitating deletion of two target genes in each tumor. Generation ofLentiviral-Cre vectors with sgRNAs targeting pairwise combinations oftumor suppressors will uncover co-operative and antagonisticinteractions between tumor suppressors in a highly-parallel manner.

Multiplexed In Vivo Genome Editing to Enable Combination TherapyScreening.

Due to the adaptive ability of many systems, combination therapies areemerging as an effective way to treat many diseases. The sheer number ofpotential therapeutic combinations quickly creates a difficult situationwhere every combination could never be tested in patients or evenpre-clinical animal models. However, this un-assayable matrix of drugcombinations could contain a combination that would work for patients.Combining drug treatments with CRISPR/Cas-mediated deletion of genescoding for additional drug targets could allow multiplexed modeling oftherapeutic combinations. Interrogation of the effects of >100 pairwisedrug treatments can be performed in parallel using the compositions andmethods described in this disclosure. For example, investigating thesepermutations in the context of three lung cancer genotypes in mousemodels of human lung cancer, would generate a semi-high-throughputsystem to interrogate the effects of pairwise drug targeting in vivo.

Extension to Other Cancer Types:

The methods described herein can be used to uncover pharmacogenomicsusceptibilities of cell growth/proliferation (e.g., in the context ofneoplasms, e.g., lung adenocarcinoma) and the methods can be applied toany convenient cancer type and/or any convenient situation in whichpopulation size of distinguishable lineages is of interest. For example,the approaches outlined in this disclosure could be adapted to anycancer that can be induced in genetically-engineered models (e.g.,sarcoma, bladder cancer, prostate cancer, ovarian cancer, pancreascancer, hematopoietic, etc.), e.g., using viral vectors.

With the wide diversity of tumor genotypes within human lungadenocarcinoma and the growing number of potential therapies, themultiplexed quantitative platform described in this disclosure canbecome a mainstay of translational cancer biology. The approachesdescribed herein will allow translational studies to effectively matchthe correct therapies to the correct patients and will have a directimpact on patient care in the clinic. It can also help with carrying outclinical trials with a subpopulation of patients that have tumors thatare the likeliest to respond to treatment—thus improving success rate ofdrug development and also rescuing drugs that had failed in a lesstargeted clinical trial.

Examples of Non-Limiting Aspects of the Disclosure

Aspects, including embodiments, of the present subject matter describedabove may be beneficial alone or in combination, with one or more otheraspects or embodiments. Without limiting the foregoing description,certain non-limiting aspects of the disclosure numbered 1-57 areprovided below. As will be apparent to those of skill in the art uponreading this disclosure, each of the individually numbered aspects maybe used or combined with any of the preceding or following individuallynumbered aspects. This is intended to provide support for all suchcombinations of aspects and is not limited to combinations of aspectsexplicitly provided below:

1. A method of measuring population size for a plurality of clonal cellpopulations in the same tissue, the method comprising:

(a) contacting a biological tissue with a plurality of cell markers thatare heritable and distinguishable from one another, to generate aplurality of distinguishable lineages of heritably marked cells withinthe contacted tissue;

(b) after sufficient time has passed for at least a portion of theheritably marked cells to undergo at least one round of division,detecting and measuring quantities of at least two of the plurality ofcell markers present in the contacted tissue, thereby generating a setof measured values; and

(c) calculating, using the set of measured values as input, a number ofheritably marked cells present in the contacted tissue for at least twoof said distinguishable lineages of heritably marked cells.

2. The method of 1, wherein the heritably marked cells within thecontacted tissue are neoplastic cells.3. The method of 1 or 2, wherein said tissue comprises neoplastic cellsand/or tumors prior to step (a).4. The method of any one of 1 to 3, wherein said detecting and measuringof step (b) is performed on a biological sample collected from thetissue.5. The method of any one of 1 to 3, wherein said detecting and measuringof step (b) is performed on a tissue sample of the contacted tissue.6. The method of any one of 1 to 5, wherein each cell marker of theplurality of cell markers corresponds to a known cell genotype for alineage of heritably marked cells.7. The method of any one of 1 to 6, wherein said contacting comprisesgenetically altering cells of the tissue to generate the heritablymarked cells.8. The method of any one of 1 to 7, wherein said method is a method ofmeasuring tumor size for a plurality of tumors of the same tissue.9. The method of any one of 1 to 8, wherein the step of contacting thetissue comprises inducing neoplastic cells.10. The method of any one of 1 to 9, wherein the cell markers are agentsthat induce or modify neoplastic cell formation and/or tumor formation.11. The method of any one of 1 to 10, wherein said detecting andmeasuring is performed after sufficient time has passed for tumors toform in the contacted tissue as a result of said contacting.12. The method of any one of 1 to 11, wherein the plurality of cellmarkers comprises barcoded nucleic acids.13. The method of 12, wherein said detecting and measuring compriseshigh-throughput sequencing and quantification of the number of sequencereads for each detected barcode.14. The method of any one of 1 to 13, wherein the plurality of cellmarkers comprises barcoded nucleic acids that induce neoplastic cellformation.15. The method of any one of 12 to 14, wherein the barcoded nucleicacids induce neoplastic cell formation and include one or more of:homology directed repair (HDR) DNA donor templates, nucleic acidsencoding one or more oncogenes, nucleic acids encoding one or morewildtype proteins, nucleic acids encoding one or more mutant proteins,nucleic acids encoding one or more CRISPR/Cas guide RNAs, nucleic acidsencoding one or more short hairpin RNAs (shRNAs), and nucleic acidsencoding one or more genome editing proteins.16. The method of 15, wherein the genome editing protein is selectedfrom: a CRISPR/Cas RNA-guided protein, a CRISPR/Cas RNA-guided proteinfused to a transcriptional activator or repressor polypeptide, a Cas9protein, a Cas9 protein fused to a transcriptional activator orrepressor polypeptide, a zinc finger nuclease (ZFN), a TALEN, aphage-derived integrase, a Cre protein, a Flp protein, and ameganuclease protein.17. The method of any one of 12 to 16, wherein the barcoded nucleicacids are linear or circular DNA molecules.18. The method of any one of 12 to 16, wherein the barcoded nucleicacids are selected from: plasmids, synthesized nucleic acid fragments,and minicircles.19. The method of any one of 12 to 16, wherein the barcoded nucleicacids are RNA molecules.20. The method of any one of 12 to 16, wherein the barcoded nucleicacids are RNA/DNA hybrids or nucleic acid/protein complexes.21. The method of any one of 1 to 19, wherein the tissue is aninvertebrate tissue.22. The method of any one of 1 to 19, wherein the tissue is a vertebratetissue.23. The method of any one of 1 to 19, wherein the tissue is a mammalianor a fish tissue.24. The method of any one of 1 to 19, wherein the tissue is a rattissue, a mouse tissue, a pig tissue, a non-human primate tissue, or ahuman tissue.25. The method of any one of 1 to 24, wherein the tissue is part of aliving animal.26. The method of any one of 1 to 24, wherein the tissue is anengineered tissue grown outside of an animal.27. The method of any one of 1 to 26, wherein the tissue is selectedfrom: muscle, lung, bronchus, pancreas, breast, liver, bile duct,gallbladder, kidney, spleen, blood, gut, brain, bone, bladder, prostate,ovary, eye, nose, tongue, mouth, pharynx, larynx, thyroid, fat,esophagus, stomach, small intestine, colon, rectum, adrenal gland, softtissue, smooth muscle, vasculature, cartilage, lymphatics, prostate,heart, skin, retina, reproductive system, and genital system.28. The method of any one of 1 to 27, wherein after sufficient time haspassed for at least a portion of the heritably marked cells to undergoat least one round of division, the method further comprises: (i)detecting and/or measuring a biomarker of the heritably marked cells,and (ii) categorizing the heritably marked cells based on the results ofsaid detecting and/or measuring of the biomarker.29. The method of 28, wherein the biomarker of one or more of: cellproliferation status, cell type, developmental cell lineage, cell death,and cellular signaling state.30. The method of any one of 1 to 29, wherein the cell markers aredelivered to the tissue via viral vector.31. The method of 30, wherein the viral vector is selected from: alentiviral vector, an adenoviral vector, an adeno-associated viral (AAV)vector, and a retroviral vector.32. A method of measuring tumor size for a plurality of clonallyindependent tumors of the same tissue, the method comprising:

(a) contacting a tissue with a plurality of barcoded nucleic acid cellmarkers, thereby generating a plurality of distinguishable lineages ofheritably marked neoplastic cells within the contacted tissue;

(b) after sufficient time has passed for at least a portion of theheritably marked neoplastic cells to undergo at least one round ofdivision, performing high-throughput nucleic acid sequencing to detectand measure quantities of at least two of the of barcoded nucleic acidcell markers present in the contacted tissue, thereby generating a setof measured values; and

(c) calculating, using the set of measured values as input, a number ofheritably marked neoplastic cells present in the contacted tissue for atleast two of said distinguishable lineages of heritably markedneoplastic cells.

33. The method of 32, wherein said tissue comprises neoplastic cellsand/or tumors prior to step (a).34. The method of 32 or 33, wherein the high-throughput nucleic acidsequencing of step (b) is performed on a biological sample collectedfrom the tissue.35. The method of 32 or 33, wherein the high-throughput nucleic acidsequencing of step (b) is performed on a tissue sample of the contactedtissue.36. The method of any one of 32 to 35, wherein each barcoded nucleicacid cell marker of the plurality of barcoded nucleic acid cell markerscorresponds to a known cell genotype for a lineage of heritably markedneoplastic cells.37. The method of any one of 32 to 36, wherein said contacting comprisesgenetically altering cells of the tissue to generate the heritablymarked neoplastic cells.38. The method of any one of 32 to 37, wherein the barcoded nucleicacids induce neoplastic cell formation.39. The method of any one of 32 to 37, wherein the barcoded nucleicacids induce neoplastic cell formation and include one or more of:homology directed repair (HDR) DNA donor templates, nucleic acidsencoding one or more oncogenes, nucleic acids encoding one or morewildtype proteins, nucleic acids encoding one or more mutant proteins,nucleic acids encoding CRISPR/Cas guide RNAs, nucleic acids encodingshort hairpin RNAs (shRNAs), and nucleic acids encoding a genome editingprotein.40. The method of 39, wherein the genome editing protein is selectedfrom: a CRISPR/Cas RNA-guided protein, a CRISPR/Cas RNA-guided proteinfused to a transcriptional activator or repressor polypeptide, a Cas9protein, a Cas9 protein fused to a transcriptional activator orrepressor polypeptide, a zinc finger nuclease (ZFN), a TALEN, aphage-derived integrase, a Cre protein, a Flp protein, and ameganuclease protein.41. The method of any one of 32 to 40, wherein the barcoded nucleicacids are linear or circular DNA molecules.42. The method of any one of 32 to 40, wherein the barcoded nucleicacids are selected from: plasmids, synthesized nucleic acid fragments,and minicircles.43. The method of any one of 32 to 42, wherein the barcoded nucleicacids are RNA/DNA hybrids or nucleic acid/protein complexes.44. The method of any one of 32 to 43, wherein the tissue is aninvertebrate tissue.45. The method of any one of 32 to 43, wherein the tissue is avertebrate tissue.46. The method of any one of 32 to 43, wherein the tissue is a mammalianor a fish tissue.47. The method of any one of 32 to 43, wherein the tissue is a rattissue, a mouse tissue, a pig tissue, a non-human primate tissue, or ahuman tissue.48. The method of any one of 32 to 47, wherein the tissue is part of aliving animal.49. The method of any one of 32 to 47, wherein the tissue is anengineered tissue grown outside of an animal.50. The method of any one of 32 to 49, wherein the tissue is selectedfrom: muscle, lung, bronchus, pancreas, breast, liver, bile duct,gallbladder, kidney, spleen, blood, gut, brain, bone, bladder, prostate,ovary, eye, nose, tongue, mouth, pharynx, larynx, thyroid, fat,esophagus, stomach, small intestine, colon, rectum, adrenal gland, softtissue, smooth muscle, vasculature, cartilage, lymphatics, prostate,heart, skin, retina, and reproductive system, and genital system.51. The method of any one of 32 to 50, wherein after sufficient time haspassed for at least a portion of the heritably marked neoplastic cellsto undergo at least one round of division, the method further comprises:(i) detecting and/or measuring a biomarker of the heritably markedneoplastic cells, and (ii) categorizing the heritably marked neoplasticcells based on the results of said detecting and/or measuring of thebiomarker.52. The method of 51, wherein the biomarker of one or more of: cellproliferation status, cell type, developmental cell lineage, cell death,and cellular signaling state.53. The method of any one of 32 to 52, wherein the cell marker isdelivered to the tissue via viral vector.54. The method of 53, wherein the viral vector is selected from: alentiviral vector, an adenoviral vector, an adeno-associated viral (AAV)vector, a bocavirus vector, a foamy virus vector, and a retroviralvector.55. The method of any one of 1-54, wherein the method includescontacting the tissue with a test compound (e.g., test drug) anddetermining whether the test compound had an effect on cell populationsize and/or distribution of cell population sizes.56. The method of any one of 1-55, wherein, after generating theheritably marked cells, the method includes transplanting one or more ofthe heritably marked cells (e.g., transplanting one or more tumors) intoone or more recipients (e.g., a secondary recipient, e.g., to seedtumors in the secondary recipient).57. The method of 56, where a test compound is administered to the oneor more recipients and the method comprises detecting and measuringquantities of at least two of the plurality of cell markers present inthe recipient(s) (e.g., to assess growth of the transplanted cells inresponse to the presence of the test compound).

EXAMPLES

The following examples are put forth so as to provide those of ordinaryskill in the art with a complete disclosure and description of how tomake and use the present invention, and are not intended to limit thescope of what the inventors regard as their invention nor are theyintended to represent that the experiments below are all or the onlyexperiments performed. Efforts have been made to ensure accuracy withrespect to numbers used (e.g. amounts, temperature, etc.) but someexperimental errors and deviations should be accounted for. Unlessindicated otherwise, parts are parts by weight, molecular weight isweight average molecular weight, temperature is in degrees Centigrade,and pressure is at or near atmospheric.

Example 1: Tuba-See: A Quantitative and Multiplexed Approach to Uncoverthe Fitness Landscape of Tumor Suppression In Vivo

Cancer growth and progression are multi-stage, stochastic evolutionaryprocesses. While cancer genome sequencing has been instrumental inidentifying the genomic alterations that occur in human tumors, theconsequences of these alterations on tumor growth within native tissuesremains largely unexplored. Genetically engineered mouse models of humancancer enable the study of tumor growth in vivo, but the lack of methodsto quantify the resulting tumor sizes in a precise and scalable mannerhas limited our ability to understand the magnitude and the mode ofaction of individual tumor suppressor genes. Here, we present a methodthat integrates tumor barcoding with ultra-deep barcode sequencing(Tuba-seq) to interrogate tumor suppressor function in mouse models ofhuman cancer. Tuba-seq uncovers different distributions of tumor sizesin three archetypal genotypes of lung tumors. By combining Tuba-seq withmultiplexed CRISPR/Cas9-mediated genome editing, we further quantifiedthe effects of eleven of the most frequently inactivatedtumor-suppressive pathways in human lung adenocarcinoma. This approachidentifies the methyltransferase Setd2 and the splicing factor Rbm10 asnovel suppressors of lung adenocarcinoma growth. With unprecedentedresolution, parallelization, and precision Tuba-seq enables a broadquantification of the fitness landscape of tumor suppressor genefunction.

Results

Tumor Barcoding with Ultra-Deep Barcode Sequencing (Tuba-Seq) Enablesthe Precise and Parallel Quantification of Tumor Sizes.

Oncogenic KRAS is a key driver of human lung adenocarcinoma, and earlystage lung tumors can be modeled using LoxP-Stop-LoxP Kras^(G12D)knock-in mice (Kras^(LSL-G12D/+)) in which expression of Cre in lungepithelial cells leads to the expression of oncogenic Kras^(G12D). LKB1and P53 are frequently mutated tumor suppressors in oncogenicKRAS-driven human lung adenocarcinomas and Lkb1- and p53-deficiencyincrease tumor burden in mouse models of oncogenic Kras^(G12D)-drivenlung tumors (FIG. 7a ). Viral-Cre-induced mouse models of lung cancerenable the simultaneous initiation of a large number of tumors andindividual tumors can be stably tagged by lentiviral-mediated DNAbarcoding. Therefore, we sought to determine whether high-throughputsequencing of the lentiviral barcode region from bulk tumor-bearinglungs could quantify the number of cancer cells within each uniquelybarcoded tumor (FIG. 7b ).

To interrogate the growth of oncogenic Kras^(G12D)-driven lung tumors aswell as the impact of Lkb1 and p53 loss on tumor growth, we initiatedlung tumors in Kras^(LSL-G12D/+); Rosa26^(LSL-Tomato) (KT), KT;Lkb1^(flox/flox) (KLT), and KT; p53^(flox/flox) (KPT) mice with alibrary of Lentiviral-Cre vectors containing greater than 10⁶ unique DNAbarcodes (Lenti-mBC/Cre; FIG. 1a and FIG. 7b ). Eleven weeks after tumorinitiation, KT mice developed widespread hyperplasias and some smalltumor masses (FIG. 1b and FIG. 7c ). Interestingly, while KLT mice hadlarge tumors of relatively uniform size, KPT mice had a very diverserange of tumor sizes (FIG. 1b ).

To quantify the cancer cell number in every lesion using ultra-deepsequencing, we PCR-amplified the integrated lentiviral barcode regionfrom ˜ 1/10^(th) of bulk lung DNA isolated from each mouse and sequencedthis to an average depth of greater than 10⁷ reads per mouse (FIG. 1a ,Methods). We observed over one-thousand-fold variation in tumor sizeswithin mice (FIG. 1c ). Barcode reads from small lesions could representunique tumors or be generated from recurrent sequencing errors ofsimilar barcodes from larger tumors. To minimize the occurrence of thesespurious tumors, we aggregated reads expected to be derived from thesame tumor barcode using an algorithm that generates a statistical modelof sequencing errors (DADA2: FIG. 2 and FIG. 8). The DADA2 aggregationrate and minimum tumor size were also selected to maximizereproducibility of our tumor-calling pipeline (FIG. 8d-f ). Theseapproaches greatly limit, but likely do not entirely eliminate, theeffect of recurrent sequencing errors on tumor quantification (FIG. 2a).

Quantification of the absolute number of cancer cells in each tumorwould allow the aggregation of data from individual mice of the samegenotype and the comparison of tumor sizes across genotypes. To enablethe conversion of read count to cancer cell number, we added cells withknown barcodes to each lung sample at a defined number prior to tissuehomogenization and DNA extraction (FIG. 1a and FIG. 9). Thus, bynormalizing tumor read counts to “benchmark” read counts we couldcalculate the absolute number of cancer cells in each tumor in eachmouse (FIG. 1a and FIG. 9).

Tuba-seq is highly reproducible between technical replicates and isinsensitive to many technical variables that could bias tumor sizedistributions including sequencing errors, variation in the intrinsicerror rate of individual Illumina® sequencing machines, barcode GCcontent, barcode diversity, tumor number within mice, and read depth(FIG. 2b-d , FIG. 10). While moderate measurement error exists at smallsizes, this does not bias the overall size distributions. Tumor sizedistributions were also highly reproducible between mice of the samegenotype (R²>0.98; FIG. 2e,f , FIG. 10g ). In fact, unsupervisedhierarchical clustering of size distributions clearly separated miceaccording to their genotype, even when tumors were induced withdifferent titers of Lenti-mBC/Cre (FIG. 2g and FIG. 10d ). Our methoddid, however, detect variation in the spectrum of tumor sizes betweenmice of the same genotypes. This variation is much greater than therandom noise observed between two fractions of tumors within the samemouse suggesting that Tuba-seq is significantly more precise than theintrinsic variability in tumor burden between mice (FIG. 2e,g ). Thus,Tuba-seq rapidly and precisely quantified the number of cancer cellswithin thousands of lung lesions in KT, KLT, and KPT mice (FIG. 1c ,FIG. 10c ).

Analysis of Tumor Sizes Uncovers Two Modes of Tumor Suppression

To assess the effect of either p53- or Lkb1-deficiency on tumor growth,we calculated the number of cancer cells in the tumors at differentpercentiles within the distribution. Interestingly, while tumors in KLTmice were consistently larger than KT tumors, deletion of p53 did notalter the number of cancer cells in the vast majority of tumors (FIG.3a-c ). Instead, a small fraction of p53-deficient tumors grew toexceptional sizes, and were among the largest in any of the mice (FIG.1c ).

To better understand the difference in tumor growth imparted by p53- andLkb1-deficiency, we defined the mathematical distributions that best fitthe tumor size distributions in KT, KLT, and KPT mice. Lkb1-deficienttumors were lognormally distributed across the full range of thedistribution (FIG. 3d ). A lognormal distribution is expected fromsimple exponential tumor growth with normally distributed rates. Toestimate average tumor size without allowing very large tumors togreatly shift this metric, we also calculated the maximum likelihoodestimator of the mean number of cancer cells given a lognormaldistribution of tumor sizes (LN mean). By this measure KLT tumors had,on average, 7-fold more cancer cells than KT tumors, consistent with therole of Lkb1 in restraining proliferation (FIG. 3a,c ). Despite greatertumor burden and visibly larger tumors in KPT mice, p53-deficiency didnot increase our estimate of mean lesion size. Instead, p53-deficienttumors were power-law distributed at large sizes and the elevated totaltumor burden was driven by rare, exceptionally large tumors (FIG. 3d ).This suggests that p53-deficient tumors acquire additional rare, yetprofoundly tumorigenic events that drive subsequent rapid growth.

Generation of a Library of Barcoded Lentiviral Vectors for MultiplexedCRISPR/Cas9-Mediated Inactivation of Tumor Suppressor Genes

Human lung adenocarcinomas have diverse genomic alterations but there isa paucity of quantitative data describing their impact on tumor growth(FIGS. 7a and 12b ). To simultaneously quantify the tumor-suppressivefunction of many known and candidate tumor suppressor genes in parallel,we combined Tuba-seq and conventional Cre-based mouse models withmultiplexed CRISPR/Cas9-mediated in vivo genome editing (FIG. 4a-c ).Assessing different tumor genotypes in a single mouse should alsomaximize the resolution of Tuba-seq, by eliminating the effect ofmouse-to-mouse variability. We first confirmed efficient Cas9-mediatedgene inactivation in lung tumors in mice with an H11^(LSL-Cas9) alleleby initiating tumors with Lentiviral-sgRNA/Cre vectors targeting eitherthe tdTomato reporter or Lkb1 (FIG. 11). Homozygous inactivation oftdTomato was achieved in around 40% of tumors and Cas9-mediated Lkb1inactivation increased tumor burden (FIG. 11). These data demonstrateour ability genetically alter tumors in Kras-driven lung cancer modelsusing these methods.

We selected eleven known and putative lung adenocarcinoma tumorsuppressor genes, which represent diverse pathways, including genes thatare broadly involved in chromatin remodeling (Setd2 and Arid1a),splicing (Rbm10), DNA damage response (Atm and p53), cell cycle control(Rb1 and Cdkn2a), nutrient and oxidative stress sensing (Lkb1 andKeap1), environmental stress responses (p53), as well as TGF-β and Wntsignaling (Smad4 and Apc, respectively) (FIG. 4b and FIG. 7a ). Weidentified efficient sgRNAs that generated indels early in thetranscripts, upstream of known functional domains, and upstream of mostmutations present in human tumors (FIG. 12a ). To allow accuratequantification of the number of cancer cells in each tumor usingTuba-seq, we diversified each tumor suppressor-targeting Lenti-sgRNA/Cre vector and four Lenti-sg Inert/Cre negative control vectors witha two-component barcode. This barcode consisted of a unique 8-nucleotide“sgID” specific to each sgRNA and a random 15-nucleotide barcode (BC) touniquely tag each tumor (sgID-BC; FIG. 4a,b and FIG. 12c-e ). In vitrocutting efficiency was determined for each of the sgRNAs individuallyand within the pool (FIG. 13).

Parallel Quantification of Tumor Suppressor Function In Vivo

To quantify the effect of inactivating each gene on lung tumor growth ina multiplexed manner, we initiated tumors in KT and KT; H11^(LSL-Cas9)(KT; Cas9) mice with a pool of the eleven barcoded Lenti-sgRNA/Crevectors and four barcoded Lenti-sg Inert/Cre vectors (Lenti-sgTS-Pool/Cre; FIG. 4b,c ). Despite receiving a lower dose of viruscompared to KT mice, KT; Cas9 mice had an increase in the number andsize of macroscopic tumors relative to KT mice 12 weeks after tumorinitiation (FIG. 4d,e ). To determine the number of cancer cells in eachtumor with each sgRNA, we amplified the sgID-BC region from bulktumor-bearing lung DNA, deep sequenced the product, and applied ourTuba-seq analysis pipeline. We calculated the entire distribution ofgrowth effects for each tumor suppressor relative to the distribution ofinert sgRNAs within each mouse. For each sgRNA, the number of cancercells in the tumors at different percentiles within the distributionwere divided by the sizes of the corresponding percentiles in the inertdistribution (FIG. 5a ). This relative and within-mouse comparisonmaximized the precision of Tuba-seq (Methods). We also determined therelative lognormal (LN) mean size of tumors containing each of theeleven tumor-suppressor-targeting sgRNAs to identify tumor suppressorsthat generally repress cancer growth (FIG. 5b ). These analysesconfirmed the known tumor-suppressive function of Lkb1, Rb1, Cdkn2a, andApc in Kras^(G12D)-driven lung tumor growth (FIG. 5a,b and FIG. 12b ).Tumors initiated with Lenti-sg TS-Pool/Cre in KT mice (which lack theH11^(LSL-Cas9) allele) had only minor differences in the sizedistributions of tumors with each sgRNA (FIG. 14a-c ).

To assess the reproducibility of this method, we analyzed an additionalcohort of KT; Cas9 mice 15 weeks after tumor initiation with Lenti-sgTS-Pool/Cre. We confirmed the tumor-suppressive effect of all tumorsuppressors identified at 12 weeks post-tumor initiation (FIG. 5c andFIG. 14e-f ). Our ability to detect tumor suppressors using multiplexedLentiviral-sgRNA/Cre delivery and tumor barcode sequencing wasreproducible as assessed by both the LN mean size and the relativenumber of cancer cells in the 95^(th) percentile tumor (FIG. 5c and FIG.14e,f ). The growth effects at the 95th percentile tumors wereexceedingly well correlated (R²=0.953) and the p-values associated withLN mean were similar between the two time points despite the use of only3 mice at the 15 week time-point (FIG. 5c ).

Identification of p53-Mediated Tumor Suppression and Recapitulation ofTumor Size Distributions within the Tumor Suppressor Pool

Consistent with the distribution of tumor sizes in KPT mice, neither LNmean nor the analysis of tumors up to the 95^(th) percentile uncoveredan effect of targeting p53 in KT; Cas9 mice with Lenti-sg TSPool/Creinitiated tumors (FIG. 5). As anticipated, Lenti-sgp53/Cre-initiatedtumors exhibited a power-law distribution at larger sizes and sgp53 wasenriched within the largest tumors in KT; Cas9 mice with Lenti-sgTSPool/Cre induced tumors (FIG. 15a,b ). This is consistent with p53inactivation enabling a small fraction of tumors to grow to large sizes.The effect of targeting p53 was greater at the later 15-week time pointconsistent with the progressive accumulation of additional alterationsand the known effect of p53 in limiting tumor progression (FIG. 15a ,FIG. 15b ).

Importantly, in KT; Cas9 mice with Lenti-sg TSPool/Cre initiated tumorsLkb1-deficient tumors exhibited a lognormal distribution of tumor sizesconsistent with the data from KLT mice (FIG. 16a ). Thus, bothp53-deficient and Lkb1-deficient tumors generated throughCRISPR/Cas9-mediated genome editing have similar size distributions tothose initiated using traditional floxed alleles. This suggests thateven in a pooled setting, quantification of individual tumor sizes canuncover distinct and characteristic distributions of tumor sizes upontumor suppressor inactivation.

Identification of Setd2 and Rbm10 as Suppressors of Lung Tumor Growth InVivo

Interestingly, in addition to appropriately uncovering several tumorsuppressors with known effects on lung tumors growth in vivo, Tuba-seqalso identified the methyltransferase Setd2 and the splicing factorRbm10 as major suppressors of lung tumor growth. Setd2 is the solehistone H3K36me3 methyltransferase and may also affect genome stabilitythrough methylation of microtubules. Despite being frequently mutated inseveral major cancer types, including lung adenocarcinoma, very littleis known about its role as a tumor suppressor in vivo. Setd2inactivation dramatically increased tumor size, with many sgSetd2-containing tumors having greater than five-fold more cancer cellsthan control tumors (FIG. 5a,b and FIG. 16b ). Interestingly, tumorsinitiated with Lenti-sg Setd2/Cre exhibited a lognormal distribution oftumor sizes (FIG. 16c ). In fact, only Lkb1-inactivation generated asimilar fitness advantage, underscoring the potential importance ofSETD2 mutations in driving rampant tumor growth in lung adenocarcinomapatients (FIG. 16).

Splicing factors have also emerged as potential tumor suppressors inmany cancer types. Although components of the spliceosome are mutated in10-15% of human lung adenocarcinomas, very little is known about theirfunctional contribution to tumor suppression. Rbm10 inactivationsignificantly increased the number of cancer cells in the top 50 percentof lung tumors and increased the LN mean size (FIG. 5a,b ). These datasuggest that the absence of Setd2-mediated lysine methylation andaberrant pre-mRNA splicing each have profound pro-tumorigenic effects inlung adenocarcinoma.

Tuba-Seq is a Precise and Sensitive Method to Quantify Tumor SuppressionIn Vivo

Quantifying the number of cancer cells in many tumors harboring distinctgenetic alterations within the same mouse allowed for the identificationand elimination of multiple sources of biological and technicalvariation (Methods). By initiating many lesions per mouse, barcodingevery lesion, pooling multiple sgRNAs into each mouse, and includinginert sgRNAs with the pool we could identify and correct for manysources of variability in tumor growth. Without these key features, ouranalysis would have been confounded by variability in the number ofinitiated tumors (CV=27%), mean tumor sizes between mice of the samegenotype (CV=38%), as well as a subtle correlation between the meaneffect size of inactivating different tumor suppressor genes withinindividual mice (CV=11%).

By calculating the size of each tumor, rather than using bulkmeasurements such as the representation of the sgRNA within all tumors,we more precisely and sensitively ascertained the growth effect ofinactivating different tumor suppressors. Interestingly, two thirds ofour identified tumor suppressors (Apc, Rb1, Rbm10, and Cdkn2a) were onlyidentified when we considered the number of cancer cells in eachbarcoded tumor, but not when we only considered the fold change in sgIDrepresentation (FIG. 5d ). In fact, effect size, statisticalsignificance, and ability to detect tumor suppressors with small effectwere all improved using the Tuba-seq pipeline compared to simplyanalyzing the change in sgID representation (FIG. 5e,f ). Thus, Tuba-seqprovides the level of resolution required to accurately capture thegrowth-suppressing effects of functional tumor suppressor genes.

Confirmation of On-Target CRISPR/Cas9-Mediated Genome Editing

As an orthogonal approach to investigate the selection for tumorsuppressor inactivation and to confirm on-target sgRNA-mediated genomeediting, we PCR-amplified and deep-sequenced each sgRNA-targeted regionfrom bulk lung DNA from three Lenti-sg TS-Pool/Cre infected (transdcued)KT; Cas9 mice. A relatively high fraction of Setd2, Lkb1, and Rb1alleles had inactivating indels at the targeted sites consistent withon-target sgRNA activity and the expansion of tumors with inactivationof these genes (FIG. 6a and FIGS. 15 c-f and 17 a,b).

Amplification and sequencing of the targeted regions of these genes frombulk lung DNA from Lenti-sg TS-Pool/Cre infected (transdcued) KT; Cas9mice also confirmed that all targeted genes contained indels (FIG. 6a ).Although all of the genes included in our pool are recurrently mutatedin human lung adenocarcinoma and frequently mutated in tumors withoncogenic KRAS (FIG. 7a ), Arid1a, Smad4, Keap1, and Atm were notidentified by any metrics as tumor suppressors (FIGS. 5 and 6 a, andFIG. 14d-f ). The lack of tumor-suppressive function of Atm isconsistent with results using an Atm^(floxed) allele, and we confirmedthe lack of tumor-suppressive function of Smad4 on oncogenicKras^(G12D)-driven lung tumor growth in vivo in KT; Cas9 mice infected(transdcued) with Lenti-sg Smad4/Cre (FIG. 17c,d ). For these genes,changes in gene expression or environmental state, additional time, orcoincident oncogene and/or tumor suppressor alterations may be requiredfor inactivation of these pathways to confer a growth advantage in lungcancer cells.

To further validate the tumor-suppressive effect of Setd2 and to assessthe histology of Setd2-deficient tumors, we induced tumors in KT and KT;Cas9 mice with lentiviral vectors containing an inert sgRNA (sgNeo2) oreither of two distinct sgRNAs targeting Setd2. KT; Cas9 mice with tumorsinitiated with either of the Lenti-sg Setd2/Cre vectors developed largeadenomas and adenocarcinoma and had significantly greater overall tumorburden than KT mice with tumors initiated with the same virus (FIG. 6b,c). While histological analysis of these mice uncovered largemouse-to-mouse variability, the analysis of individual tumor sizes byTuba-seq confirmed a nearly four-fold increase in the number of cancercells in Setd2-deficient tumors relative to control tumors (FIG. 6c,dand FIG. 18). Importantly, the validation of Setd2 tumor suppression byconventional methods required many more mice than our initial screen ofeleven putative tumor suppressors emphasizing the benefit ofmultiplexing sgRNAs to increase throughput and decrease costs.

Discussion

While many putative tumor suppressors have been identified from cancergenome sequencing, limited strategies exist to test their function(e.g., in vivo) in a rapid, systematic, and quantitative manner (FIG.19). By combining DNA barcoding, high-throughput sequencing, andCRISPR/Cas9-mediated genome editing, Tuba-seq not only increases thethroughput of these analyses, but also enables exceptionally precise anddetailed quantification of tumor growth in vivo.

Interestingly, tumors initiated at the same time, within the same mouse,with the same genomic alterations grew to vastly different sizes afteronly 12 weeks of growth. Thus, additional spontaneous alterations,differences in the state of the initial transformed cell, or the localmicroenvironment may impact how rapidly a tumor grows and whether it hasthe capacity for continued expansion. Tuba-seq was also uniquely able touncover genotype-specific distributions of tumor sizes that revealedproperties of gene function. p53-deficiency generated a tumor sizedistribution that is power-law distributed for the largest tumors,consistent with a Markov process where very large tumors are generatedby additional, rarely acquired driver mutations. Conversely, Lkb1inactivation increased the size of a majority of lesions suggesting anordinary exponential growth process. Thus, tumor suppressors can havedifferent modes of tumor suppression, identified via Tuba-seq, that mayportend their molecular function. Interestingly, Setd2 has recently beensuggested to methylate tubulin, and Setd2-deficiency can lead to variousforms of genomic instability including micronuclei and laggingchromosomes due to alterations in microtubules. Genome instability wouldbe expected to generate rare, advantageous alterations and tumors growththat is highly-stochastic and power-law distributed. However, the sizedistribution of Setd2-deficient lung tumors in our studies was strictlylognormal, therefore we speculate that the main impact of Setd2 loss isthe induction of gene expression programs that generally dysregulategrowth (FIG. 6d and FIG. 16b,c ).

The scale of our analyses, which evaluated thousands of individualtumors, dramatically improved our ability to identify functional tumorsuppressor genes. Estimating tumor growth via bulk measurements wouldhave identified only a third of the tumor suppressors that we uncoveredas advantageous to tumor growth (FIG. 5d-f ). Unlike conventional floxedalleles, CRISPR/Cas9-mediated genome editing in the lung generatedhomozygous null alleles in approximately half of all tumors (FIG. 11d ).Thus, while the lack of uniform homozygous deletion of targeted geneswould reduce the tumor suppressive signal from bulk measurements, bybarcoding and analyzing each tumor, Tuba-seq effectively overcomes thistechnological limitation.

By analyzing a large number of tumor suppressors, our data suggest thatearly neoplastic cells reside in an evolutionarily nascent state wheremany tumor suppressor alterations were adaptive and conferred a growthadvantage. In contrast, CRISPR/Cas9 screens in cancer cell lines havefound that additional tumor suppressor alterations provide littleadvantage and can even be detrimental. This finding is consistent withcancer cell lines residing in a much more mature evolutionary state,approaching optimal growth fitness due to their origin fromadvanced-stage disease as well as selection for optimal proliferativeability in culture. Furthermore, the intimate link between tumorsuppression and many aspects of the in vivo environment underscores theimportance of analyzing the effects of tumor suppressor loss in tumorsin vivo (or for example in the context of a tissue such as an organoidculture or 3D cultured tissue).

Interestingly, the frequency of tumor suppressor alterations in humancancer did not directly correspond to the magnitude of their tumorsuppressor function. For example, SETD2 and RBM10 are mutated in similarpercentages of human lung adenocarcinomas, but Setd2-deficiencyconferred a much greater growth benefit than Rbm10-deficiency (FIG. 5a,b). This highlights the growing need for methods that allow rapid andquantitative analyses of gene function in vivo to determine thefunctional importance of low-frequency putative tumor suppressors thatmay be profoundly important for individual patients.

There is a very limited understanding of the clinical importance oftumor suppressor alterations, and this remains a major unmet need, butstrong drivers of tumor growth may represent more attractive clinicaltargets than weak drivers. Tuba-seq permits investigation of morecomplex combinations of tumor suppressor gene loss, as well as theanalysis of other aspects of tumor growth and progression. Tuba-seq isalso adaptable to study other cancer types and should allow theinvestigation of genes that normally promote, rather than inhibit, tumorgrowth. Finally, this method allows the investigation ofgenotype-specific therapeutic responses which could ultimately lead tomore precise and personalized patient treatment.

Statistical Properties of Lesions in this Study

The distributions of tumor sizes were generally lognormal withinclinations towards a 2^(nd)-order power law when looking within aMouse-sgRNA pair (FIG. 20). Each tumor in our study was assigned alog-transformed size t_(mrb) defined by the mouse m that harbored it,the cognate sgRNA r identified by its first barcode, and a uniquebarcode sequence (consensus of the DADA2 cluster) b. Our approach wasdesigned to interrogate and address a variety of sources of error: wefound that (i) the number of instigated tumors within replicate mice,(often littermates) infected (transdcued) with the same lentiviral titervia the same intubation procedure, varied greatly, (ii) the mean tumorsize varied within replicate mice, (iii) certain mice were more amenableto growth of tumors with sgRNAs targeting specific tumor suppressors,and (iv) the size of tumors with the same sgRNA within the same mousevaried dramatically.

CV Interrogated by Minimized using Source of Variance Efficiency ofInfection 27% Random DADA2 (transduction) Barcodes clustering Mouse 38%Inert sgRNAs Normalization to Inert Mean Mouse-Tumor Suppressor 11%Multiplexed PCA Mixture Interactions sgRNAs Model Stochastic Progression511%  Multiple Lognormal MLE Infections of Mean Variable of InterestTumor Suppressor 31% n.a.

Overall, the effect of a tumor suppressor inactivation on tumor burdenis small compared to these other sources. Prior viral-Cre-basedgenetically engineered mouse models address the main source ofvariability by initiating hundreds to thousands of tumors per mouse. Weobserve stochasticity in the size of tumors instigated within the samemouse with the same genetic constructs even in this setting. The numberof cancer cells in individual tumors in these experiments is nevermeasured accurately; instead, total tumor area is most often measured,which is a conflation of mean tumor size and the number of instigatedtumors. Thus, this approach is flawed because (i) the sampled mean sizeis not the best estimator of mean size, (ii) the number of instigatedtumors is never directly measured (a quantity that varies with aCoefficient of Variation (CV) of 27%), (iii) variability in the mousebackground is ignored, and (iv) the methods used to assess tumor areaalso introduce variability. For these reasons, the magnitude of effectof even the most powerful tumor suppressor in a Kras^(G12D/+) background(Setd2) is smaller than the variance between replicate mice (FIG. 6c ).

The variability in the number of lesions instigated by viral-Cre vectorswithin individual mice also affects estimates of a tumor suppressor'seffect. By uniquely barcoding each tumor and then precisely callingtumors using our computational approach detailed in the Methods, weminimize this source of variability. Variance in the number of calledlesions in our pipeline exhibited a CV of 10.7% across repeatedsequencing runs, whereas the variance in called lesions betweenreplicate mice exhibited a CV of 27%. Therefore, our estimates of tumornumber, based on unique DNA barcodes, are significantly more precisethan presuming the number of tumors is constant between replicate mice(which likely have different numbers of epithelial cells infected(transdcued) by lentiviral vectors due to technical variability with theinfection (transduction)). Below, we interrogate, mitigate, and discussthe remaining sources of variability listed in the above Table.

sgRNA-Agnostic Mouse-to-Mouse Variability

Our multiplexed approach interrogates growth effects attributable to (i)the CRISPR/Cas9-target tumor suppressor gene, (ii) the individual mice,and (iii) their interactions. This is only possible because we includedmany sgRNAs within each mouse and measured many lesions with the samesgRNA in the same mouse. We observed a statistically significantdifference between ostensibly replicate mice in the mean,log-transformed, bias-corrected expectation size of each sgRNA(η_(mr)=E_(mr)[t_(mrb)]) These differences could be succinctlysummarized and then subtracted from t_(mrb) to better resolve thestrength of each tumor suppressor.

Mice exhibited both sgRNA-agnostic growth perturbationsη_(m)=E_(r)[η_(mr)] (i.e. there was a spectrum of tumor-susceptible andtumor-resistive mice) and sgRNA-dependent covariance within the miceη_(mr) (e.g. mice that harbored unusually large Lkb1-deficient tumorsalso harbored unusually large Setd2-deficient tumors). About 40% of themouse-to-mouse variability was eliminated by correctly normalizingη_(m), while the variability in η_(mr) not ascribable to sgRNA-agnosticfactors was estimated to be only 10.7%. We could only eliminate a fifthof this η_(mr) variability (detailed below). Therefore, most of thevariability in tumor susceptibility appears to be sgRNA-agnostic,however subtle gene-mouse covariance is still consequential whenestimating average tumor growth advantages to precisions <10%.

Replicate mice, i.e. those with the same genetic-engineered elementsanalysed at the same time-point after tumor initiation, were oftenlittermates and cage-mates, but descend from a mixed 129/BL6backgrounds. While these mice likely have a far more homogenous genotypeand environment than real-world patients, relevant differences betweenindividual mice still emerged. It is important to note that while thesetrends can be identified in our data due to our unprecedentedresolution, the variation is small and should have an even greatereffect on experiments that compare different mice constructs (forexample conventional approaches that compare tumor growth in mice withand without a floxed allele of a gene of interest, or our own resultsfrom mice with tumors initiated with Lenti-sg Setd2/Cre versusLenti-sgNeo/Cre (see FIG. 6)).

Because each mouse contained several inert sgRNAs (whose means did notdiffer appreciably from each other within a mouse) we were able tosubtract the sgRNA-agnostic mice effects simply by normalizing sizesrelative to the aggregated inert sgRNA mean:μ_(mr)=E_(mr)[t_(mrb)]−E_(m,inerts)[t_(mrb)]. In our nonparametricapproach, we simply divide by the median inert sgRNA, which tends to bealmost identical the LN MLE mean.

sgRNA-Specific Mouse-to-Mouse Variability

The availability of multiple active sgRNAs in a single mouse allowed usto interrogate sgRNA-specific mouse effects. Overall, the μ_(mr) matrixwas highly positively correlated between active sgRNAs in KT; Cas9 mice.We decomposed these correlations using Principle Component Analysis(PCA). The first Principle Component (PC1) explained 75% of the variancein μ_(mr) across replicate KT; Cas9 mice. We tested several hypothesesfor this covariance:

-   -   1. Mice that harbored larger tumors on average might also harbor        larger tumor variance in log-scale. If so, then sg Lkb1 to sg        Inert tumor size ratio would co-vary with the sg Setd2 to sg        Inert tumors size ratio.    -   2. Mouse gender drives these diverging growth patterns.    -   3. Cas9 endonuclease cutting efficiency varied between mice that        are H11^(LSL-Cas9/+) versus mice that are        H11^(LSL-Cas9/LSL-Cas9).    -   4. An unknown genetic or environmental factor within the mouse        perturbs the strength of a subset of drivers.

We investigated these first two hypotheses by comparing PC1 to meantumor size and by comparing KT; H11^(LSL-cas9/+) to KT;H11^(LSL-Cas9/LSL-Cas9) mice in our KT; Cas9 12-week cohort. PC1correlated well with both mean tumor size (as calculated via ourpipeline) and lung weight (FIG. 20b-d ). Lung weight (in grams) wasdetermined at the time of lung sample collection and is likelyinfluenced by tumor number and mean tumor size. The correlation of lungweight with PC1, like mean tumor size, ensures that these observedtrends are not a pipeline artefact. Mouse gender also co-varied with PC1(Point-Biserial Correlation r=0.75, data not shown) and is consistentwith our first hypothesis, as male mice exhibit both larger tumors and alarger size discrepancy between strong drivers and inerts.

H11^(LSL-Cas9) allele status (heterozygous or homozygous) was notstatistically-significantly correlated with PC1 (r=0.34, data not shown)in 12-week KT; Cas9 mice. Therefore, we do not believe that beingheterozygous or homozygous for the H11^(LSL-Cas9) allele contributessubstantially the efficacy of gene inactivation.

Finally, the hypothesis of a latent genetic or environmental factor istoo open-ended to be tested here. However, our methodology permitsstudies of these factors moving forward.

Thus, we conclude that tumor permissivity and mouse gender are mostlyresponsible for these sgRNA-specific differences between replicate mice,and that Cas9 endonuclease heterozygosity does not seem to appreciablyaffect tumor growth, and results from our analysis pipeline concur withother mouse measurements.

A Mixture of Probabilistic Principle Components model was used toeliminate μ_(mr) from μ_(mr). This model defines the log-likelihood of amouse arising from the same distribution as the others in its cohort ofreplicates. In essence, this model identifies mice with anomalous sgRNAprofiles. However, rather than categorize mice as either ‘outlier’ or‘acceptable’ mice, we simply weighted each mouse based its likelihood ofoutlying. Statistically, an ‘outlier’ is defined as a point that appearsto be drawn from a different distribution than its cohort. Indeed, wefound that similar outlier mice were identified using Mahalanobisdistance—a common metric for identifying outliers in multidimensionaldata. However, the Mahalanobis distance metric requires some thresholdfor classifying outliers that would be ad hoc in our application.Weighting mice using our Mixture of Probabilistic Principle ComponentsModel, reduced the variability of E_(r)[μ_(mr)] for KT; Cas9 mice by2.1%. Although this is only a mild improvement, we included thiscorrection in our final report of the mean growth advantage conferred byan sgRNA because we felt that this value should account for every sourceof variability identified. The final reported mean growth effect of eachsgRNA in a cohort of replicate mice was an arithmetic mean of μ_(mr)across all mice weighted by the likelihood of each mouse m in ourmixture model P(m; μ_(mr)), i.e.

$\frac{\Sigma_{m}\mspace{14mu} {P\left( {m;\mu_{mr}} \right)}\mu_{mr}}{\Sigma_{m}\mspace{14mu} {P\left( {m;\mu_{mr}} \right)}}.$

Our Parametric and Nonparametric Approaches, and Statistical Tests

Comprehensively measuring the size spectrums of tumor growth and thenidentifying all the exogenous factors behind this spectrum presents aconundrum: growth advantages could be summarized with a highly-processedmeasure of tumor size that accounted for every known, quantifiedconcern, or growth advantages could be summarized more explicitly in amanner that makes fewer assumptions. We chose both extremes. Thequalitative conclusions do not differ much in either case; however, wepresent both approaches because the agreement is encouraging and becausethe different approaches may appeal to readers with differentsensibilities.

Our approach based on Maximum Likelihood estimation is detailed in thesection above. Summarily, it attempts to account for our understandingof the lognormal shape of size distributions, and mouse-to-mousevariability in (i) the number of instigated tumors, (ii) overall tumorpermissivity, and (iii) sgRNA-specific variability. It leverages themulti-dimensionality of our size measurements and corrects for everyknown exogenous factor that we found. Below, we discuss the limitationsof assuming log-normality and extend the parametric approach to tumorsuppressor distributions exhibiting power-law tails.

Our nonparametric summary presents percentiles of t^((nonparametric))_(mrb) (defined above) to assay for increased tumor growth at variouslocations in the size distributions. It makes no assumption of the shapeof the tumor distributions and does not model mouse-to-mousevariability. Although by correcting for the median size of inerts andthe number of tumors residing in each mouse, a majority of themouse-to-mouse variability is eliminated. For this reason, after thefirst experiment, percentiles were always reported relative to theircorresponding inert percentile. Autocorrelation between differentpercentile tiers for the sgRNA is expected and observed; the differentpercentile tiers are not statistically-independent values and we deployno statistical test that assumes their independence.

The nonparametric approach generally finds that the 90 to 99^(th)percentiles of distributions of active sgRNAs are maximally deviant fromthe inerts. Our finding that distributions are at least log-normallyskewed is consistent with this phenomenon. Furthermore, active sgRNAscan introduce inframe insertions and deletions that should mimic theinert distribution, so we expected the smallest tumors in an activesgRNA distribution—with in-frame mutations or no mutations—to mimicinert sizes. Lastly, the haploinsufficiency of a single null allele isgenerally unknown, but if haploinsufficiency is partially dominant ornon-existent then size distributions would be most deviant at higher (90to 99^(th)) percentiles.

Therefore, we used the 95^(th) percentile as a crude summary of thegrowth benefit of a driver, as it approximately balanced our concerns ofthe null-mutation rate, zygosity, statistical resolution (which declinesat higher percentiles), and our understanding of the size distributionsin general. Our data suggest that loss of a tumor suppressor does notnecessarily lead to a growth advantage across all individual tumors (forexamples p53-versus Lkb1-deficiency in FIGS. 1 and 2). Indeed, the 95thpercentile measure fails to detect p53 in our experiment for reasonsthat are in line with the expected consequence of p53 loss andfat-tailed distributions. Nonetheless, simplifications can be useful andthe 95th percentile of sizes summarizes differences in growth well.

All confidence intervals and p-values were obtained via bootstrapping oft_(mrb). After bootstrap sampling, all subsequent steps in our analysispipelines were recalculated for every bootstrap (normalizations toinerts, PCA, etc.). Bootstrapping samples were equal in size to theoriginal t_(mrb) for each experiment (e.g. the tumors in KT; Cas9 miceanalyzed 12-weeks after tumor initiation) and were sampled withreplacement. 200,000 samples were drawn for every 95% confidenceinterval reported and 2,000,000 samples were drawn for every p-valuereported. Confidence intervals of ratios reflect uncertainty in both theactive sgRNA distribution and the inert sgRNA distribution. Therefore,when the confidence interval of an sgRNA ratio does not subsume 1, thenull hypothesis that this summary statistic of the sgRNA matches theinert sgRNA can be rejected with p<0.05 (assuming no correction forMultiple Hypotheses).

All p-values report the two-sided hypothesis that an sgRNA summarystatistic differs from the inert sgRNA summary statistic and wereBonferroni-corrected for our multiple hypotheses that any one of the 11active sgRNAs could incur a growth advantage or disadvantage. Whileactive sgRNAs were always compared to the entire sgRNA distribution(four different inert sgRNAs), the inert sgRNAs were compared only tothe distribution of the other three inerts. p-values were not reportedbeyond 0.0001, as this is the resolution limit of bootstrapping whenlimited to 2,000,000 samples.

Comprehensive Parametric Description of Size Distributions

Lesion sizes were approximately lognormally distributed with excessivequantities of very large lesions in some genotypes. We fit a widevariety of 2-3 parameter probability distributions to the observeddistribution of lesion sizes for each genotype and time: (Log)-Normal,(Log)-Gamma, (Log)-Logistic, Exponential, Beta, Generalized-ExtremeValue (including Gumbel), and Power-Law (including Pareto). All lesionsize distributions were best fit with either a Lognormal, Log-gamma, orLog-Logistic distribution, although no single distribution outperformedall others. A Kolmogorov-Smirnov test often rejected the best-fittingsingle distribution—i.e. only a least-improper fit could be found inmany cases. This shortcoming, underscores the enormous quantities oftumor sizes that we were able to measure for the first time and thecomplexities of tumor progression. Therefore, we investigatedmulti-family parametric fits.

A combination of Lognormal and Power-Law scaling, for somedistributions, best described our data. Although Log-Gamma andLog-Logistic fits were sometimes superior to Lognormal fits, thesealternative distributions merely have faster-growing higher moments,which is suggestive of Power-Law behaviour. Moreover, the MaximumLikelihood Estimators of Log-Gamma and Log-Logistic distributionparameters must be solved numerically without guarantee of convergence.

Care was taken to identify Power-Law distributions impartially.Potential Power-Law fits were parameterized using maximum likelihood andadjudicated using marginal likelihood:

-   -   1. The Maximum Likelihood Lognormal fit

${Max}_{\mu_{r},\sigma_{r}}\left\lbrack {\sum\limits_{m,b}{{\mathcal{L}}\left( {{t_{mrb};\mu_{r}},\sigma_{r}} \right)}} \right\rbrack$

for each sgRNA distribution for the entire support of positive realnumbers was identified. Here,

denotes the probability density of a lognormal distribution.

-   -   2. The Maximum Likelihood Power-Law fit for tumors over the        domain [x^((min)) _(r), ∞) was identified:

${Max}_{x_{r}^{(\min)},\alpha_{r}}\left\lbrack {\sum\limits_{m,b}{\begin{matrix}{{\mathcal{L}}\left( {{t_{mrb};\mu_{r}},\sigma_{r}} \right)} \\{{\mathcal{L}}\left( {{t_{mrb};x_{r}^{(\min)}},\alpha_{r}} \right)}\end{matrix}\mspace{14mu} {if}\mspace{14mu} \begin{matrix}{t_{mrb} < x_{r}^{(\min)}} \\{t_{mrb} \geq x_{r}^{(\min)}}\end{matrix}}} \right\rbrack$

Here

denotes the probability density of a Power-Law or Pareto Distributionwith exponent α_(r) and the lognormal fit from step 1 is used. Note thata Power-Law is undefined when x^((min))=0 and so it is customary to testpower-laws over a limited support with a freely-floating minimum.

-   -   3. The multi-fit model was adjudicated using Marginal        Likelihood: the likelihood of the observed data corrected for        the model's degrees of freedom using Bayesian-Information        Criterion (BIC).

This approach is recommended by Alstot et al and their accompanyingsoftware package was used for this analysis. Details of the maximumlikelihood fits are provided in FIG. 3. Reported p-values are atransformation of the Marginal Likelihoods of a joint lognormal andpower-law fit, such that p=1/(1+Exp[Marginal Likelihood]). These valuestest the null hypothesis that the data is lognormal-distributedthroughout its entire support.

We also test the hypothesis that sizes are distributed according to anExponentially-Truncated Power Law. This comparison is a commoncounter-hypothesis to an ordinary Power Law and would imply thatscale-free behavior exists only over a finite interval³. We do not seegood evidence for Exponentially-truncated Power Law behaviour (data notshown). For this reason, we believe that the data support scale-freemodels of tumor progression in the genotypes discussed below.

We observe strong, recurring evidence that p53-deficient tumors arepower-law distributed at large scales. Power-law dynamics were observedin all incarnations of the Kras^(G12D/+/p)53Δ genotype (KPT tumors, andboth KT; Cas9 sgp53 tumor time-points). The marginal likelihoods for allof these Power-Law distribution fits were good or excellent. Thisagreement strongly supports the hypothesis that the Kras^(G12D/+)/p53Δgenotype is power-law distributed in tumor size.

In general, the ML exponent of power-law fits was approximately two(α˜2). Power-law dynamics have been hypothesized to explain cancerincidence rates, however tumor size distributions have not been wellstudied because measurements to test this hypothesis were previouslyprohibitively time-consuming and costly, so we explored a simpleevolutionary model that yields a power-law distribution of sizes in thenext section. Our deep interrogation of lesion sizes proved useful innot only precisely identifying driver growth advantages, but also inuncovering aspects of their underlying mode of action.

Additional Rarely-Acquired Driver Mutations Predict a Power LawDistribution of Tumor Sizes

p53-deficient tumors exhibit a Power-Law distribution of sizes in theirrightmost tail (FIG. 3d ). Power law distributions generally do notarise from a single-step Markov process and, instead, arise fromcompound random processes, e.g. random walks or accretion processes⁶.The simplest, and we believe most-likely, explanation for this observedpower law distribution is a combination of exponential processes, namelythe rare acquisition of a second driver event inexponentially-expanding, p53-deficient tumors.

Suppose that tumor size N(t) initially grows exponentially over time twith rate r₁, such that N(t)=e^(r) ¹ ^(t). Let N(t=0)=1, i.e. there isone tumorigenic cell at infection (transduction) time which is definedas t=0. Furthermore, we assume that at time t* a new clone with a newdriver emerges in the tumor population and grows at a much faster rater₂, such that this clone dominates the tumor population at the time ofsacrifice t^(F), i.e. r₂(t^(F)−t*)>>r₁t^(F). Note that 0≤t*<t^(F).Lastly, suppose that this transformative clone emerges randomly in timewith a probability that is proportional to the size of the tumor, i.e.p(t*)˜μN(t). In this scenario, the size of tumors at time of analysisN(t=t^(F))=n is

n=e ^(r2(t) ^(F) ^(−t*))

n∝e ^(−r) ² ^(t*)

Based on the derivation reviewed in Newman, M. Power laws, Paretodistributions and Zipf's law. Contemp. Phys. 46, 323-351 (2005),entitled Combinations of exponentials (Section 4.1), we find:

${p(n)} = \frac{n^{- {({1 + \frac{r_{1}}{r_{2}}})}}}{r_{2}}$

Tumor sizes are power-law distributed with exponent

$1 + {\frac{r_{1}}{r_{2}}.}$

This result implies either that the observed exponent must be less than2 or that additional drivers must be acquired. The Maximum-Likelihoodestimate of the exponent for KPT mice sacrificed at 11 weeks is slightlygreater than two, while the exponent for sgp53 tumors in KT; Cas9 micesacrificed at 15 weeks is slightly less than two (although both of thesevalues subsume two in their 95% CI).

Assumption Description N(t) = e^(rt) Exponential growth dynamicsp(t*)~μN(t) 2^(nd) driver arises w/probability proportional topopulation size r₂(t^(F) − t*) >> r₁t^(F) 2^(nd) driver completesselective sweep by time of sacrifice

All of the above assumptions made in other basic mathematical models oftumor progression⁴. Thus, we believe that a Markov Process is the bestexplanation of the observed Power-Law.

Lastly, we note that the transformative event at time (is unspecified.It could be a genetic alteration, an epigenetic change, a switch in cellsignalling state, etc. We further note that there are other processesthat may generate a Power-Law distribution.

Evidence of Tumors with Multiple Lentiviral Infections

Size measurements are precise enough to identify lesions putativelyinfected (transdcued) by multiple lentiviral vector. Our firstexperiment (KT, KLT, KPT m ice) used larger viral titters (6,000 to22,000 capsids), so we expected multiple infections to be more common.If two different viral vectors infected (transdcued) the same foundingcell, then it would expand into a single tumor annotated as twolesions—by both lentiviral barcodes. Therefore, if we observed twobarcoded tumors of the same size within an individual mouse, then wemight expect that these arose from two lentiviral vectors initiating asingle lesion. Thus, we investigated the size difference between eachlesion and its nearest neighbor in the same mouse.

We observed that a small fraction of lesions were closer in size thanexpected, suggesting that some lesions may have arisen from cells thatwhere initiated by infection (transduction) with more than onelentiviral vector. Our (null) expected distribution represents the sizedifferences between an observed lesion size and their nearest neighborin different (randomly-selected) mice. Although our data suggests thatmultiple infections occur in ˜1% of tumors, we do not believe that thatthis rare occurrence substantially affects the other conclusions of ourstudy because (i) multiple infections appear to be rare, and (ii)multiple infections should attenuate our estimates of a driver's growthbenefit (as multiple-infections would confer a growth advantage to ourbaseline—the sg Inert constructs). Nonetheless, this preliminarydiscovery once again illustrates our approaches ability to uncover newbiology using old techniques.

Methods

Mice and Tumor Initiation

Kras^(LSL-G12D) (K), Lkb1^(flox) (L), p53^(flox) (P), R26^(LSL-Tomato)(T),and H11^(LSL-Cas9) (Cas9) mice have been described. Lung tumors wereinitiated by intratracheal infection (transduction) of mice aspreviously described using lentiviral-Cre vectors at the titersindicated. Tumor burden was assessed by fluorescence microscopy, lungweight, and histology as indicated. All experiments were performed inaccordance with Stanford University Institutional Animal Care and UseCommittee guidelines.

Generation of Barcoded Lenti-mBC/Cre and Lenti-sg Pool/Cre Vector PoolsTo enable quantification of the number of cancer cells in individualtumors in parallel using high-throughput sequencing, we diversifiedlentiviral-Cre vectors with a short barcode sequence that would beunique to each tumor by virtue of stable integration of the lentiviralvector into the initial infected (transdcued) lung epithelial cell. Wegenerated tumors in a variety of mouse backgrounds with two differentpools of barcoded lentiviral vectors. The first was a pool of ˜10⁶uniquely barcoded variants of Lenti-PGK-Cre (Lenti-millionBC/Cre;Lenti-mBC/Cre, generated by pooling six barcoded Lenti-U6-sgRNA/PGK-Crevectors) which we used to analyze the number of cancer cells in tumorsinduced in Kras^(LSL-G12D/+); R26^(LSL-Tomato) (KT), Kras^(LSL-G12D/+);p53^(flox/flox); R26^(LSL-Tomato) (KPT), andKras^(LSL-G12D/+; Lkb)1^(flox/flox); R26^(LSL-Tomato) (KLT) mice (FIG.1). The second was a pool of 15 barcoded Lenti-U6-sgRNA/PGK-Cre vectorswhich we used to assess the tumor suppressive effect of candidate tumorsuppressor genes in three different genetic backgrounds by infecting KT;H11^(LSL-Cas9) (KT; Cas9) and KT mice. Our Lenti-sg Inert/Cre vectorsincluded three sgRNAs that target the NeoR gene within theRosa26^(LSL-Tomato) allele, which are actively cutting, but functionallyinert, negative control sgRNAs.

Design, Generation, and Screening of sgRNAs

We generated lentiviral vectors carrying Cre as well as an sgRNAtargeting each of 11 known and putative lung adenocarcinoma tumorsuppressors: sg Lkb1, sgP53, sgApc, sgAtm, sgArid1a, sgCdkn2a, sgKeap1,sg Rb1, sgRbm10, sg Setd2, and sgSmad4. Vectors were also generatedcarrying inert guides: sgNeo1, sgNeo2, sgNeo3, sgNT1, and sgNT3. Allpossible 20-bp sgRNAs (using an NGG PAM) targeting each tumor suppressorgene of interest were identified and scored for predicted on-targetcutting efficiency using an available sgRNA design/scoring algorithm¹⁰.For each tumor suppressor gene, we selected three unique sgRNAspredicted to be the most likely to produce null alleles; preference wasgiven to sgRNAs with the highest predicted cutting efficiencies, as wellas those targeting exons conserved in all known splice isoforms(ENSEMBL), closest to splice acceptor/splice donor sites, positionedearliest in the gene coding region, occurring upstream of annotatedfunctional domains (InterPro; UniProt), and occurring upstream of knownhuman lung adenocarcinoma mutation sites. Lenti-U6-sgRNA/Cre vectorscontaining each sgRNA were generated as previously described. Briefly,Q5 site-directed mutagenesis (NEB E0554S) was used to insert sgRNAs intothe parental lentiviral vector containing the U6 promoter as well asPGK-Cre. The cutting efficiency of each sgRNA was determined byinfecting LSL-YFP; Cas9 cells with each Lenti-sgRNA/Cre virus.Forty-eight hours after infection (transduction), flow cytometricquantification of YFP-positive cells was used to determine percentinfection (transduction). DNA was then extracted from all cells and thetargeted tumor suppressor gene locus was amplified by PCR.

PCR amplicons were Sanger sequenced and analyzed using TIDE analysis toquantify percent indel formation. Finally, the indel percent determinedby TIDE was divided by the percent infection (transduction) of LSL-YFP;Cas9 cells, as determined by flow cytometry, to determine sgRNA cuttingefficiency. The most efficient sgRNA targeting each tumor suppressorgene of interest was used for subsequent experiments. sgRNAs targetingTomato and Lkb1 have been described previously, and we previouslyvalidated an sgRNA targeting p53 (unpublished data). Primers sequencesused to amplify target indel regions for the top guides used in thisstudy are below:

F primer R primer (5′→3′) (5′→3′) sgApc _1 TGACTTTGCA CCCACTCCCCGGGCAAGTTT TGTTACCTTT sgAridla_3 CAGCAGTCCC GGAGCCATTT CAACTCCATACTTGGGGTTA sgAtm_3 GCCCCAAGTG AGCTCTGGCT AGAATCAGTG CCTTGTGGATsgCdkn2a_2 GGCTTCTTTC GGCTCATTTG TTGGGTCCTG GGTTGCTTCT sgKeap1 _2CTGAGCCAGC GGCCTATCCC AACTCTGTGA ACTTCTGAGC sgRb1 _3 AACTGTGCTGACACCACCAC GTGTGTGCAA CACCATCATC sgRbm10_3 CAAAGCTGGA CTGGCTGGAGAGCGAGACTG CTGTGAGAGT sgSetd2_1 TCTGCAAGTT TGGATTCAGGT CAAGCGATGAAGCCTAGATGG sgSetd2_2 CCTCCAGCCG GAACGCCGAA CTCCTCAT CCTAAGCAG sgSmad4_3GCCTTTCTGT TTCCAGGCTG GGAAATGGAA AGTGGTAAGG sgNeo_1 TTGTCAAGACCCACCATGAT CGACCTGTCC ATTCGGCAAG sgNeo_2 TCTGGACGAA GCTCCAATCCGAGCATCAGG TTCCATTCAA sgNeo_3 CGCTGTTCTC TGGATACTTTC CTCTTCCTCATCGGCAGGA

Barcode Diversification of Lenti-sgRNA/Cre

After identifying the best sgRNA targeting each tumor suppressor ofinterest, we diversified the corresponding Lenti-sg RNA/Cre vector witha known 8-nucleotide ID specific to each individual sgRNA (sgID; bold)and the 15-nucleotide random barcode (BC; underlined) (see FIG. 10a ).

Primer (5′→3′) Universal AGCTAGGGATCCGC Reverse CGCATAACCAGTG PrimerBarcoded AGCTAGTCCGGNNNNN Forward NNNAANNNNNTTNNNN PrimerNAANNNNNATGCCCAA GAAGAAGAGGAAGGTGTC

These primers were used to PCR amplify a region of the Lenti-PGK-Crevector that included the 3′ end of the PGK promoter and the 5′ part ofCre. PCR was performed using PrimeSTAR® HS DNA Polymerase (premix)(Clontech, R040A) and PCR products were purified using the Qiagen® PCRPurification Kit (28106). The PCR insert was digested with BspEI andBamHI and ligated with the Lenti-sgRNA-Cre vectors cut with XmaI (whichproduces a BspEI compatible end) and BamHI.

To generate a large number of uniquely barcoded vectors, we ligated 300ng of each XmaI, BamHI-digested Lenti-sgRNA-Cre vector with 180 ng ofeach BspEI, BamHI-digested PCR product using T4 Ligase (NEB, M0202L) andstandard protocols (80 μl total reaction volume). Ligations were PCRpurified using the Qiagen® PCR Purification Kit to remove residual salt.To obtain a pool of the greatest possible number of uniquely barcodedLenti-sgRNA/Cre vectors, 1 μl of purified ligation was transformed into20 μl of ElectroMAX DH10B cells (Thermo Fisher, 18290015). Cells wereelectroporated in 0.1 cm GenePulser/MicroPulser Cuvettes (Bio-Rad,165-2089) in a BD MicroPulser™ Electroporator (Bio-Rad, 165-2100) at 1.9kV. Cells were then rescued by adding 500 μl media and shaking at 200rpm for 30 minutes at 37° C. For each ligation, bacteria were plated onseven LB-Amp plates (1 plate with 1 μl, 1 plate with 10 μl, and 5 plateswith 100 μl). The following day, colonies were counted on the 1 μl or 10μl plate to estimate the number of colonies on the 100 μl plates, andthis was used as an initial estimation of number of unique barcodesassociated with each ID.

10 ml of liquid LB-Amp was added to each plate of bacteria to pool thecolonies. Colonies were scraped off of the plates into the liquid, andall plates from each transformation were combined into a flask. Flaskswere shaken at 200 rpm for 30 minutes at 37° C. to mix. DNA wasMidi-prepped using the Qiagen® HiSpeed MidiPrep Kit (12643). DNAconcentrations were determined using a Qubit dsDNA HS Kit (Invitrogen,Q32851).

As a quality control measure, the sgID-BC region from eachLenti-sgRNA-sgID-BC/Cre plasmid pool was PCR amplified with GoTaq Greenpolymerase (Promega M7123) following manufacturer's instructions. ThesePCR products were Sanger sequenced (Stanford PAN facility) to confirmthe expected sgID and the presence of a random BC. Since BspEI and XmaIhave compatible overhangs but different recognition sites, theLenti-sgRNA-sgID-BC/Cre vectors generated from successful ligation ofthe sgID/BC lack an XmaI site. Thus for pools that had a detectableamount of unbarcoded parental Lenti-sgRNA/Cre plasmid as determined bySanger sequencing (>5%), we destroyed the parental unbarcoded vector bydigesting the pool with XmaI (NEB, 1000 reaction) using standardmethods. These re-digested plasmid pools were re-purified using theQiagen® PCR Purification Kit and concentration was redetermined byNanoDrop.

Generation of Lenti-mBC/Cre and Lenti-TS-Pool/Cre To obtain a librarywith approximately 10⁶ associated barcodes to use in our initialexperiments in mice that lacked the H11^(LSL-Cas9) allele, we pooled sixsgID-BC barcoded vectors to create Lenti-million Barcode/Cre(Lenti-mBC/Cre). We then pooled the barcoded Lenti-sgRNA-sgID-BC/Crevectors (sgLkb1, sgp53, sgApc, sgAtm, sgArid1a, sgCdkn2a, sgKeap1,sgNeo1, sgNeo2, sgNeo3, sgNT1, sgRb1, sgRbm10, sgSetd2, and sgSmad4) togenerate Lenti-sg TS-Poot/Cre. All plasmids were pooled at equal ratiosas determined by Qubit concentration prior to lentivirus production.

Production, Purification, and Titering of Lentivirus

Lentiviral vectors were produced using polyethylenimine (PEI)-basedtransfection of 293T cells with the lentiviral vectors and delta8.2 andVSV-G packaging plasm ids. Lenti-mBC/Cre, Lenti-sgTS-Pool/Cre, Lenti-sgTomato/Cre, Lenti-sg Lkb1, Lenti-sg Setd2#1/Cre, Lenti-sg Setd2#3/Cre,Lenti-sgNeo2/Cre, and Lenti-sg Smad4/Cre were generated for tumorinitiation. Sodium butyrate (Sigma Aldrich, B5887) was added at a finalconcentration of 0.2 mM eight hours after transfection to increaseproduction of viral particles. Virus-containing media was collected 36,48, and 60 hours after transfection, concentrated by ultracentrifugation(25,000 rpm for 1.5-2 hours), resuspended overnight in PBS, and frozenat −80° C. Concentrated lentiviral particles were titered by infectingLSL-YFP cells (a gift from Dr. Alejandro Sweet-Cordero), determining thepercent YFP-positive cells by flow cytometry, and comparing theinfectious titer to a lentiviral preparation of known titer.

Generation of “Benchmark” Cell Lines

Three uniquely barcoded Lenti-Cre vectors with the sgID “TTCTGCCT” wereused to generate benchmark cell lines that could be spiked into eachbulk lung sample at a known cell number to enable the calculation ofcancer cell number within each tumor. Plasmid DNA from individualbacterial colonies was isolated using the Qiagen® QIAprep Spin MiniprepKit (27106). Clones were Sanger sequenced, lentivirus was produced asdescribed above, and LSL-YFP cells were infected (transdcued) at a verylow multiplicity of infection (transduction) such that approximately 3%of cells were YFP-positive after 48 hours. Infected (transdcued) cellswere expanded and sorted using a BD Aria II™ (BD Biosciences).YFP-positive sorted cells were replated and expanded to obtain a largenumber of cells. After expansion, cells were re-analyzed for percentYFP-positive cells on a BD LSR II™ analyzer (BD Biosciences). Using thispercentage, the number of total cells needed to contain 5×10⁵ integratedbarcoded lentiviral vectors was calculated for each of the three celllines and cells were aliquoted and frozen based on this calculation.

Summary of all Mouse Infections

Genotype Virus Type Viral Titer KT Lenti-mBC/Cre 6.8 × 10⁵ KT_(low)Lenti-mBC/Cre 1.7 × 10⁵ KPT Lenti-mBC/Cre 1.7 × 10⁵ KLT Lenti-mBC/Cre1.7 × 10⁴ KT Lenti-TS-Pool/Cre 9.0 × 10⁴ KT; Cas9 Lenti-TS-Pool/Cre 2.2× 10⁴ KT; Cas9 Lenti-sgNeo2/Cre  9 × 10³ KT; Cas9 Lenti-sgSetd2#1/Cre  9× 10³ KT; Cas9 Lenti-sgSetd2#2/Cre  9 × 10³ KT Lenti-sgSmad4/Cre 10⁵ KT;Cas9 Lenti-sgSmad4/Cre 10⁵

Isolation of Genomic DNA from Mouse Lungs

For experiments in which barcode sequencing was used to quantify thenumber of cancer cells in each tumor the whole lungs from each mousewere homogenized using a Fisher TissueMeiser. 5×10⁵ cells from each ofthe three individually barcoded benchmark cell lines were added at thetime of homogenization. Tissue was homogenized in 20 ml lysis buffer(100 mM NaCl, 20 mM Tris, 10 mM EDTA, 0.5% SDS) with 200 μl of 20 mg/mlProteinase K (Life Technologies, AM2544). Homogenized tissue wasincubated at 55° C. overnight. To maintain accurate representation ofall tumors, DNA was phenol-chloroform extracted and ethanol precipitatedfrom ˜ 1/10^(th) of the total lung lysate using standard protocols. Forlungs weighing less than 0.3 grams, DNA was extracted from ˜⅕^(th) ofthe total lung lysate, and for those weighing less than 0.2 grams, DNAwas extracted from ˜ 3/10^(th) of the total lung lysate to increase DNAyield.

Preparation of sgID-BC Libraries for Sequencing

Libraries were prepared by amplifying the sgID-BC region from 32 μg ofgenomic DNA per mouse. The sgID-BC region of the integratedLenti-sgRNA-BC/Cre vectors was PCR amplified using one of 24 primerpairs that contain TruSeq Illumina® adapters and a 5′ multiplexing tag(TruSeq i7 index region indicated in bold):

Primer (5′→3′) Universal AATGATACGGCGACCACCGAGATCTA Forward PrimerCACTCTTTCCCTACACGACGCTCTTC CGATCTGCGCACGTCTGCCGCGCTG Reverse PrimerCAAGCAGAAGACGGCATACGAGATN NNNNNGTGACTGGACTTCAGACGTGTGCTCTTCCGATCCAGGTTCTTGCG AACCTCAT

We used a single-step PCR amplification of sgID-BC regions, which wefound to be a highly reproducible and quantitative method to determinethe number of cancer cells in each tumor. We performed eight 100 μl PCRreactions per mouse (4 μg DNA per reaction) using OneTaq 2× Master Mixwith Standard buffer (NEB, M0482L) with the following PCR program:

1. 94 C 10 min

2. 94 C 30 sec

3. 550 30 sec

4. 680 30 sec

5. GO TO 2 (34×)

6. 68 C 7 min

7. 4 C infinity

Pooled PCR products were isolated by gel electrophoresis and gelextracted using the Qiagen® MinElute Gel Extraction kit. Theconcentration of purified PCR products from individual mice wasdetermined by Bioanalyzer (Agilent Technologies) and pooled at equalratios. Samples were sequenced on an Illumina® HiSeq to generate 100 bpsingle-end reads (ELIM Biopharmaceuticals, Inc).

Identifying Distinct sgRNAs and Tumors Via Ultra-Deep Sequencing

The unique sgID-BC identifies tumors. These sgID-BCs were detected vianext generation sequencing on the Illumina® HiSeq. The size of eachtumor, with respect to cell number, was expected to roughly correspondto the abundance of each unique sgID-BC pair. Because tumor sizes variedby factors larger than the read sequencing error rate, distinguishingtrue tumors from recurrent read errors required careful analysis of thedeep-sequencing data.

Tumors and their respective sgRNAs were identified in three steps: (i)abnormal and low quality reads were discarded from the ultra-deepsequencing runs, (ii) unique barcode pileups were bundled into groupsthat we predicted to arise from the same tumor, and (iii) cell numberwas estimated from these bundles in the manner that proved mostreproducible.

Read Pre-Processing

Reads contained a two-component DNA barcode (an 8-nucleotide sgID and a21-nucleotide barcode sequence that contains 15 random nucleotides) thatbegan 49 nucleotides downstream of our forward primer and ended 22nucleotides upstream of the end of our 100-bp single-end reads. Wediscarded unusual reads: those that lacked the flanking lentiviralsequences, those that contained unexpected barcodes, and those with higherror rates. This was accomplished in three steps (FIG. 8a ):

-   -   1. We examined the 12 lentiviral nucleotides immediately        upstream and downstream of the sgID-BC. These 12 nucleotides        were identified using pairs of adjacent 6-mer search strings,        such that each 6-mer could tolerate one mismatch. Although we        expected these 12 nucleotides to begin at position 37 within the        read, we did not require this positioning or leverage this        information. A nested 6-mer approach (with two opportunities to        identify the lentiviral sequences flanking the sgID-BC) was used        to minimize read discarding. This was particularly important in        this first step because the non-barcoded regions of our reads        were used to estimate sequencing error rates and, therefore,        should not be biased against read errors. For ˜7-8% of reads,        this 2^(nd) 6-mer match salvaged the read, i.e. the 6-mers        immediately flanking the sgID-BC were not as expected (despite        our tolerance of one mismatch) yet the 6-mers immediately        outside of these inner 6-mer sequences were recognizable and        allowed us to salvage the read and identify the barcodes.    -   2. We then discarded reads in which the sgID-BC deviated in        length by greater than two nucleotides in either direction.        Because our first barcode was expected to contain one of the 15        sgIDs, we discarded reads that did not match one of these 15        sequences. One mismatch and one indel were permitted in the        matching.    -   3. We then end-trimmed each read such that 18 bp flanked either        end of the sgID-BC. We then filtered the trimmed reads according        to quality score, retaining those that were predicted to contain        no more than two sequencing errors. We also discarded reads with        uncalled bases in the second (random) barcode and rectified        uncalled bases elsewhere.

In these three stages, 14% of reads were discarded at stage one, ˜7% atstage two, and <2% at stage three.

We then examined those reads that failed at each stage. By performingBLAST searches, we determined that those reads discarded at stage oneoften contained uninformative sequences corresponding to artifacts fromeither our preparation (Phi X bacteriophage genome and mouse genome) orother samples paired with us on the lane (common plasmid DNAs). In stagetwo, we found that reads with aberrant barcode lengths often containedlarge indels or had one or both of their sgID-BC completely missing.Lastly, very few reads were discarded in stage three due to the factthat internal regions of the reads exhibited higher quality scores thanthe corresponding termini. As a consequence of this trend, it is commonpractice to end-trim reads prior to discarding those reads predicted tocontain greater than two sequencing errors.

Clustering of Unique Read Pileups Via DADA2

sgID-BC reads were aggregated into sets of identical sequences andcounted. The counts of unique DNA barcode pairs do not directlycorrespond to unique tumors because large tumors are expected togenerate recurrent sequencing errors (FIG. 8b ). We therefore spentconsiderable effort developing a method to distinguish small tumors fromrecurrent sequencing errors arising from large tumors (consider, forexample, that a tumor of 10 million cells will produce sequencing-errorpileups that mimic a 10-100 thousand-cell tumor, if the Illumina®machine has a 0.1-1% error rate). DADA2 has been used previously toaddress this issue in barcoding experiments involving ultra-deepsequencing. However, because it was designed for ultra-deep sequencingof full-length IIlumina amplicons, we had to tailor and calibrate it forour purposes.

In DADA2, the likelihood of barcode pileups resulting from a recurrentsequencing error of a larger pileup depends upon:

-   -   1. The abundance of the larger pileup,    -   2. The specific nucleotide differences between the smaller and        larger pileups, and    -   3. The average quality scores of the smaller pileup at the        variant positions.

Factors one and two are, at first, considered heuristically (to maximizecomputational speed) and then more precisely (when needed) via aNeedleman-Wunsch algorithm. DADA2 splits a cluster into two when theprobability that a smaller pileup was generated by sequencing errors isless than a Therefore, this value represents a threshold for splittinglarger clusters. When this threshold is large, read pileups are splitpermissively (many called tumors, perhaps dividing large tumors), andwhen n is small, read pileups are split restrictively (few calledtumors, perhaps aggregating distinct small tumors).

The likelihood of sequencing errors was inferred from our ultra-deepsequencing data. Phred quality scores provide a theoretical estimate ofsequencing error rates, however these estimates tend to vary fromIllumina® machine to Illumina® machine and do not account for thespecifics of our protocol (including, for example, occasional errorsintroduced via PCR amplification, despite our use of high-fidelitypolymerase). Ordinarily, DADA2 will estimate sequencing error ratessimultaneously with the unique DNA clusters; however, our lentiviralconstructs had non-degenerate regions outside of our sgID-BC region thatcould be used to estimate sequencing error rates directly. Moreover,estimating error rates and barcode clusters jointly was morecomputationally intensive, requiring greater than 20,000 CPU-hours forclustering our entire dataset and exploring the relevant clusteringparameters. A sequencing error model was trained to each Illumina®machine by:

-   -   1. Generating training pseudo-reads by concatenating the 18        nucleotides immediately upstream of our sgID-BC with the 18        nucleotides immediately downstream of the barcodes, then    -   2. Clustered these pseudo-reads using a single run of DADA2.    -   3. Using the error rates estimated from this training run to        cluster the sgID-BC using a single run of DADA2.

We used a very low value of Ω=10⁻¹⁰⁰ to estimate sequencing errors inthe training run, as we expected only one cluster of lentiviralsgID-BC-flanking sequences. Altering this value did not affect ourtraining run appreciably, but we nonetheless observed occasional verysmall derivative clusters from our lentiviral sequence even at thisvalue. These derivative clusters are presumably rare DNA artifacts andnever amounted to >2% of our processed reads. We felt that using a verystringent DADA2 run to estimate sequencing errors represented a superiorapproach (by virtue of the Goldilocks principle): a more permissivethreshold might over-fit sequencing errors and underestimate sequencingerror rates, while an approach where error rates were estimated directlyfrom each read's deviance from expectation (akin to a DADA2 run wheren=0) would ignore the presence DNA artifacts in our data and, therefore,overestimate sequencing error rates.

We trained sequencing error rates on each Illumina® machine used in thisstudy (seven in total). Training allowed the probability of everysubstitution type (A-C, AST, etc) to be estimated. The error rates as afunction of Phred quality score were determined using LOESS regressionof the available data (FIG. 8c ). In general, error rates wereapproximately two to three times higher than predicted by the Phredquality scores for transversions (and approximately consistent withexpectations for transitions). This elevated error rate is typical andmay reflect miscalibration of the machines and/or mutations introducedduring PCR.

We then clustered the dual barcodes that passed our pre-processingfilters using DADA2. Barcodes were given seven nucleotides ofnon-degenerate lentiviral flanking regions so that any indels within thebarcodes could be identified (without adequate flanking sequences, DNAalignment algorithms sometimes miscall indels as multiple pointmutations). During clustering, we also required (i) that clustersdeviate from each other by at least two bases (i.e.MIN_HAMMING_DISTANCE=2), (ii) that new clusters only be formed whenpileup size exceeded expectations under the error process by at least afactor of two (MIN_FOLD=2), and (iii) that the Needleman-Wunschalgorithm consider only alignments with at most four net insertions ordeletions (BAND_SIZE=4, VECTORIZED_ALIGNMENT=FALSE). None of thesechoices affected the results appreciably, but they increasedcomputational performance and offered additional verification thatbarcodes were aggregated into tumors of reasonable size.

Vetting and Calibration of Pipeline

We sequenced our first PCR-amplified, multiplexed DNA libraries (fromKT, KLT, and KPT tumors) in triplicate to vet and design ourtumor-calling approach.

Reproducibility was interrogated in three ways: (i) the correlationbetween estimated cell abundances for all barcodes and all mice, (ii)the variation in the number of lesions called for each sgID in eachmouse in our first experiment, and (iii) the variation in mean size foreach sgID—which should be constant in mice not expressing Cas9. Becausethe read depth of our triplicate run naturally varied (40.1×10⁶,22.2×10⁶, and 34.9×10⁶ reads after pre-processing), these three runswere performed on distinct Illumina® machines with different sequencingerror rates, and, because our initial lentiviral pool contained sixdifferent sgIDs with varying levels of barcode diversity, the technicalvariability in our vetting process well-approximated the technicalvariability of later experiments. In our tumor-size analysis pipeline,we found:

-   -   1. The mean abundance of our three “benchmark” DNA barcodes was        more reproducible between replicate runs than the median        abundance. Thus, this mean value of benchmark read abundance        (corresponding to 500,000 cells) was used to convert read        abundance into the absolute cell number of cancer cells in each        tumor (FIG. 9).    -   2. Ignoring reads with ≥2 errors from the consensus barcode of a        cluster improved reproducibility. Typically, ˜80-90% of reads in        a barcode cluster were exact matches to the consensus barcode,        while ˜5% of reads were single errors from this read, and ˜5-15%        of reads deviated at ≥2 errors. These reads with ≥2 errors were        poorly correlated between replicate runs and hampered our        ability to reproducibly estimate absolute cell number/tumor        size. These reads, presumably, have neither enough evidence to        be considered their own lesion, nor sufficient evidence to be        counted towards the larger cluster. Therefore, these reads were        excluded.    -   3. The cluster-splitting proclivity of DADA2 was thresholded at        Ω=10⁻¹⁰ and required that lesions contain 500 cells for FIGS.        1-3 and ≥1000 cells for FIGS. 4-6 to maximize reproducibility        between replicate runs (FIG. 8d-f ). Threshold parameters with        high specificity (small Ω, high minimum cell number) called        lesion sizes more reproducibly, whereas threshold parameters        with high sensitivity (large Ω, low minimum cell number) called        lesion quantities more reproducibly. Over-prioritizing only one        facet of reproducibility would be imprudent. With two        thresholds, considering different facets of measurement error,        we better balanced these competing priorities.

With this pipeline, we interrogated the diversity of the barcode in ourscreen in several ways. First, we confirmed that nucleotides in thisbarcode were evenly distributed among A's, T's, C's, and G's (FIG. 10b). Second, we found no evidence for an excess of repeated string (e.g.sequences AAAAA). Third, we calculated the number of random barcodespaired to each sgID in our lentiviral pool. Due to the large number ofuniquely barcoded variants of each vector that we generated through ourbarcode ligation approach, (see Barcode diversification ofLenti-sgRNA/Cre) most barcodes that exist in our lentiviral pool werenever detected in any lesions in any of the experiments (becausediversity is much higher that total lesion number). Nonetheless, westill inferred the amount of barcode diversity from the observedbarcodes. To make this inference, we assumed that the probability ofobserving a barcode in i mice is Poisson distributed: P(k=i; λ)=λ^(k)e^(−λ)/k!, where λ_(r)=L_(r)/D_(r), is a ratio of the number of calledlesions L_(r) for each sgID r in our entire dataset (a known quantity)divided by the total number of unique barcodes D_(r) for each sgID. Bynoting that λ_(r)/(1−e^(−λr))μ_(non-zero), where μ_(non-zero)=Σi=1 ^(∞)P(k=i; λ_(r)) is simply the mean number of occurrences of each barcodethat occurred once or more, we could calculate Dr. Across our entiredataset, the average probability of the same barcode initiating twodistinct tumors in the same mouse was 0.91%.

Good barcode diversity is also demonstrated by the six sgIDs inLenti-mBC/Cre experiment. If barcode diversity was low and barcodesoverlapped often within a mouse, then the mean sizes of the less diversesgIDs would increase—as two distinct tumors with the same barcode wouldbe bundled together. However, the mean sizes of each sgID vary by <1%within replicate mice, thus refuting this possibility. We also assessedour ability to call sgIDs accurately, despite sequencing errors, byprocessing deep-sequencing runs in two ways: by identifying each read'scognate sgID before clustering based on the raw read sequence or byidentifying cognate sgIDs after clustering based on the consensussequence of the cluster. Using either approach, 99.8% of reads paired tothe same cognate sgID, thus providing assurance that sgIDs areaccurately identified. We opted to employ the latter approach for ourfinal analysis.

By thoroughly developing and vetting our tumor-calling pipeline, wesalvaged an extra decade of size resolution. Our three DNA benchmarks(added to the lung samples at the very beginning of DNA preparation)(FIG. 9) offer a glimpse of this resolution. Sequencing errors of theDNA benchmarks are easily identified by the DNA benchmark's unique sgIDand known secondary barcodes. While these sequencing errors are usuallydiscarded, we can treat them as ordinary read pileups and observe theproperties of potential sequencing errors. Without our calibratedanalysis pipeline, the sequencing errors appear as lesions of ˜10³cells; with our pipeline, these sequencing errors emerge as lesions of˜10² cells—below our minimum cell threshold (FIG. 2a ).

More importantly, our pipeline is robust to technical perturbations. Wemore intensively profiled reproducibility with two additional technicalperturbations in two specific mice from the first experiment. First, aKLT 11-week mouse (JB1349) was sequenced at great depth and thenrandomly down-sampled ten-fold to typical read depth (this down-samplingwas more dramatic than any variability in read depth actually detectedthroughout our study). Lesion sizes were very highly correlated in thisfirst perturbation (FIG. 2b ). Additionally, a KT 11-week mouse (IW1301)was amplified in two PCR reactions with different multiplexing tags(FIG. 2c ). PCR and multiplexing appears to hamper reproducibility morethan read depth, although reproducibility is good overall. These micealso display two encouraging reproducibility trends: (i) largerlesions/tumors were most consistent between replicates, and (ii) theoverall shape (histogram) of tumor lesion sizes were better correlatedbetween the replicates than individual tumors (e.g. r=0.89 for eachlesion in IW1301, whereas r=0.993 for the abundance of tumors within the60 histogram bins of FIG. 8b ). This second observation implies that ourtechnical perturbations introduce unbiased noise. Also, all correlationscompare logarithmic size; because larger tumors are better correlated,this transformation substantially reduces the Pearson correlationcoefficient.

Minimizing the Influence of GC Amplification Bias on Tumor-Size CallingWe define each tumor in our study by a size T_(mrb) corresponding to themouse m that harbored it, the cognate sgRNA r identified by its firstbarcode, and a unique barcode sequence (consensus of the DADA2 cluster)b. Given the approximately lognormal structure of our data (FIG. 3d anddata not shown), we log-transformed and normalized sizes such thatT_(mrb)=Ln(T_(mrb)/E_(mr)[T_(mrb)]). Here E_(mr)[T_(mrb)]=Σ_(b)T_(mrb)/N_(mr) is the expected lesion size for a given mouse m and sgRNArand we will use this notation for expectation values. Thisnotation—where aggregated indices are dropped from subscripts—is usedthroughout. GC biases were subtle: the coefficient of variation (CV) ofE_(mr)[T_(mrb)] was 5.0%. This marginal distribution still exhibited asubtle dependence on the GC-content of the combined barcode sequencethat was best described by a 4^(th)-order least-squares polynomial fitf₄ (b) of E_(b)[τ_(mrb)] (adjusted r²=0.994). The sgIDs were alldesigned with well-balanced GC-content, however the second barcodecomprises random sequences. While the multinomial process of generatingbarcodes made intermediate levels of GC-content most common, somedeviation of GC-content was observed. Maximal values of f₄ (b) arise atintermediate GC-content, suggesting that PCR biases amplificationtowards template DNA of intermediate melting temperature. We subtractedthe effects of this GC-bias from log-transformed values:t_(mrb)=Ln[T_(mrb)]−f₄ (b). This correction alters tumor sizes by 5% onaverage.

Calculation of In Vitro Cutting Efficiency Using the Lenti-TS-Pool/CreVirus

Cas9 expressing cell lines were infected (transdcued) withLenti-TS-Pool/Cre virus and harvested after 48 hours. gDNA was extractedand targeted loci were amplified using the above primers.

Analysis of Indels at Target Sites

To confirm CRISPR/Cas9-induced indel formation in vivo, the targetedregion of each gene of interest was PCR-amplified from genomic DNAextracted from bulk lung samples using GoTaq Green polymerase (PromegaM7123) and primer pairs that yield short amplicons amenable topaired-end sequencing:

 F primer R primer (5′→3′) (5′→3′) Apc CATGGCATAAAG TCTCCTGAACGCAGTTACTACA GCTGGATAC Aridla CCAGTCCAATG GGGTACCCATG GATCAGATG TCCTTGTTGAtm CACCCAGTTGA CCGTTTTCGGA CCCTATCTTC AGTTGACAG Cdkn2a CAACGTTCACGACCAGCGTGTC TAGCAGCTC CAGGAAG Keap1 GGCTTATTGAG GCTGCTGCACG TTCGCCTACAAGGAAGT Rb1 GGTACCCGATC AAGGAACACAG ATGTCAGAGA CTCCCACAC Rbm10TACTCAGCCGC GAGGATTTGTT TTTCTTTGC CCGCATCAG Setd2 CTGTTGTGGTTTTTTCAGTTTG GTGCCAAAG AGAACAGCCTTT Smad4 TCGATTCAAAC CTTGTGGAAGCCATCCAACA CACAGGAAT Lkb1 GGGCCTGTACC TGTCCCTTGCT CATTTGAG GTCCTAACA p53CATCACCTCAC CAGGGGTCTCG TGCATGGAC GTGACAG Neo1 GGCAGGATCTC AGTACGTGCTCCTGTCATCT GCTCGATG Neo2 CGGACCGCTAT GAGCGGCGATA CAGGACATA CCGTAAAG Neo3GATCGGCCATT CATCAGAGCAG GAACAAGAT CCGATTGT

PCR products were either gel-extracted or purified directly using theQiagen® MinElute kit. DNA concentration was determined using the QubitHS assay, following manufacturer's instructions. All 14 purified PCRproducts were combined in equal proportions for each mouse. TruSeqIIlumina® sequencing adapters were ligated on to the pooled PCR productswith a single multiplexing tag per mouse using SPRIworks (BeckmanCoulter, A88267) with standard protocols. Sequencing was performed onthe Illumina HiSeq to generate single-end, 150-bp reads (StanfordFunctional Genomics Facility).

Custom Python scripts were used to analyze the indel sequencing data.For each of the 14 targeted regions, an 8-mer was selected on eitherside of the targeted region to generate a 46 basepair region. Reads wererequired to contain both anchors and no sequencing errors were allowed.The length of each fragment between the two anchors was then determinedand compared to the expected length. Indels were categorized accordingto the number of basepairs inserted or deleted.

The percent of indels for each individual locus in each individual mousewas calculated as follows:

${\% {Indels}} = \frac{{{Total}\mspace{14mu} {Reads}} - {{WildType}\mspace{14mu} {Reads}}}{{Total}\mspace{14mu} {Reads}}$

Then the average % of indels in the three Neo loci was calculated andthe % indels at every other targeted locus was normalized to this valueto generate the % Indels relative to Neo that are plotted in FIG. 6 a.

Calculation of In Vitro Cutting Efficiency Using the Lenti-TS-Pool/CreVirus

Cas9 expressing cell lines were infected (transdcued) withLenti-TS-Pool/Cre virus and harvested after 48 hours. gDNA was extractedand targeted loci were amplified using the above primers (see Analysisof indels at target sites). First, all primers were pooled and 15 roundsof PCR were performed using GoTaq Green polymerase (Promega M7123).These products were then used for subsequent amplification withindividual primer pairs as described above. Sequencing libraries wereprepared as described above.

Histology, Immunohistochemistry, and Tumor Analysis

Samples were fixed in 4% formalin and paraffin-embedded.Immunohistochemistry was performed on 4 μm sections with the ABCVectastain kits (Vector Laboratories) with antibodies against Tomato(Rockland Immunochemicals, 600-401-379), Smad4 (AbCam, AB40759) and Sox9(EMD Milipore, AB5535). Sections were developed with DAB andcounterstained with haematoxylin. Haematoxylin and eosin staining wasperformed using standard methods.

Sections from lungs infected (transdcued) with Lenti-sg Tomato/Cre, werestained for Tomato and tumors were scored as positive (>95% Tomatopositive cancer cells), Negative (no Tomato-positive cancer cells), ormixed (all other tumors). Tumors were classified and counted from asingle section through all lung lobes from 4 independent mice.

Quantification of Tumor Area and Barcode Sequencing of Tumors Inducedwith Lenti-Sg Setd2 and Lenti-sgNeo

Tumor-bearing lung lobes from mice infected (transdcued) with Lenti-sgSetd2#1/Cre, Lenti-sg Setd2#2/Cre or Lenti-sgNeo2/Cre virus wereembedded in paraffin, sectioned, and stained with haematoxylin andeosin. Percent tumor area was determined using ImageJ.

The distribution of the number of cancer cells in individual tumors inKT; Cas9 mice infected (transdcued) with Lenti-sg Setd2#1/Cre andLenti-sgNeo2/Cre was assessed by Illumina® sequencing of theirrespective lentiviral barcodes and subsequent analysis as describedabove.

Western Blotting for Lkb1 and Cas9 Microdissected Tomato-positive lungtumors from KT and KT; Cas9 mice with Lenti-sg Lkb1/Cre initiated tumorswere analyzed for Cas9 and Lkb1 protein expression. Samples were lysedin RIPA buffer and boiled with LDS loading dye. Denatured samples wererun on a 4%-12% Bis-Tris gel (NuPage) and transferred onto a PVDFmembrane. Membranes were immunoblotted using primary antibodies againstHsp90 (BD Transduction Laboratories, 610419), Lkb1 (Cell Signaling,13031P), Cas9 (Novus Biologicals, NBP2-36440), and secondaryHRP-conjugated anti-mouse (Santa Cruz Biotechnology, sc-2005) andanti-rabbit (Santa Cruz Biotechnology, sc-2004) antibodies.

Survival Analysis of Mice with Cas9 Mediated Inactivation of Smad4

To confirm lack of functional tumor suppression attributable to Smad4,KT and KT; Cas9 mice were infected (transdcued) intratracheally with 10⁵Lenti-sg Smad4/Cre. Mice were sacrificed when they displayed visiblesigns of distress to assess survival.

Example 2: Multiplexed Quantitative Analysis of Oncogenic Variants InVivo

Large-scale genomic analyses of human cancers have catalogued somaticpoint mutations thought to initiate tumor development and sustain cancergrowth. However, determining the functional significance of specificalterations remains a major bottleneck in our understanding of thegenetic determinants of cancer. Here, we present a platform thatintegrates multiplexed AAV/Cas9-mediated homology-directed repair (HDR)with DNA barcoding and high-throughput sequencing to simultaneouslyinvestigate multiple genomic alterations in de novo cancers in mice.Using this approach, we introduced a barcoded library of non-synonymousmutations into hotspot codons 12 and 13 of Kras in adult somatic cellsto initiate tumors in the lung, pancreas, and muscle. High-throughputsequencing of barcoded Kras^(HDR) alleles from bulk lung and pancreasuncovered surprising diversity in Kras variant oncogenicity. Rapid,cost-effective, and quantitative approaches to simultaneouslyinvestigate the function of precise genomic alterations in vivo willuncover novel biological and clinically actionable insights intocarcinogenesis.

Results

To analyze the oncogenic function of diverse point mutations in vivo ina quantitative and relatively high-throughput manner, we developed aplatform for somatic AAV/Cas9-mediated HDR that incorporates DNAbarcoding and high-throughput sequencing in autochthonous mouse modelsof several cancer types (FIG. 23a-d ). We designed, generated, andvalidated a library of AAV vectors to introduce all possible Kras codon12 and 13 single-nucleotide non-synonymous point mutations into somaticmouse cells in a multiplexed manner (FIG. 23e-g and FIG. 27). Each AAVcontained an sgRNA targeting the second exon of Kras, a ˜2 kb Kras HDRtemplate, and Cre-recombinase (AAV-Kras^(HDR)/sg Kras/Cre; FIG. 23e andFIG. 27a-c ).

The Kras^(HDR) template contained either wild type (WT) Kras or one ofthe 12 single-nucleotide non-synonymous mutations in codons 12 and 13 ofKras, as well as the genomic sequence flanking the second exon of Kras.Each Kras^(HDR) template also contained silent mutations within the sgKras target sequence and associated protospacer adjacent motif (PAM*) toprevent Cas9-mediated cleavage of Kras^(HDR) alleles. To enable theparallel quantification of individual tumors by high-throughputsequencing of DNA from bulk tissue, we diversified each Kras^(HDR)template with a random eight-nucleotide barcode engineered into thewobble positions of the codons downstream of 12 and 13 (FIG. 23e andFIG. 27b,c ).

The AAV vectors also encoded Cre-recombinase. Cre-expression enabledtumor initiation in mice containing a Cre-regulated Cas9 allele(H11^(LSL-Cas9)), a fluorescent Cre-reporter allele (R26^(LSL-Tomato)),as well as floxed alleles of the well-known tumor suppressor genes p53(p53^(flox)) or Lkb1 (Lkb1^(flox)). We packaged the AAV-Kras^(HDR)/sgKras/Cre library using an AAV8 capsid that enables high titerproduction, efficient transduction of mouse lung epithelial cells invivo (FIG. 28), and transduction of a wide range of adult mousetissues³⁵.

We initially transduced Cas9-expressing cells in culture withAAV-Kras^(HDR)/sg Kras/Cre to determine whether AAV/Cas9-mediated HDRwould be an unbiased method to engineer point mutations into theendogenous Kras locus (FIG. 27e ). Kras^(HDR)-specific PCR amplificationfollowed by high-throughput sequencing of transduced cells confirmed thegeneration of all point mutant Kras alleles (FIG. 27f,g ). Furthermore,in vitro Kras^(HDR) allele frequencies correlated with theirrepresentation in the AAV-Kras^(HDR)/sg Kras/Cre plasm id library. Thisresult confirms that HDR using our AAV vector is not discernably biasedby any single-nucleotide Kras codon 12 or 13 point mutation in theKras^(HDR) template. Therefore, any differential expansion of tumorsharboring specific Kras mutant alleles can be attributed to biochemicaldifferences between Kras variants rather than differences in theefficiency of HDR using donor DNA templates with each Kras allele (FIG.27h ).

To determine whether HDR in somatic cells could initiate tumors, and toinvestigate whether Kras variants differ in their ability to drivetumorigenesis, we delivered AAV-Kras^(HDR)/sg Kras/Cre libraryintratracheally to the lungs of mice with the H11^(LSL-Cas9) allele(FIG. 24 and FIG. 29). Specifically, we transduced three differentgenotypes of mice to provide insight into whether concurrentinactivation of tumor suppressor genes modulates Kras variantoncogenicity: 1) Rosa26^(LSL-Tomato); H11^(LSL-Cas9) (T; H11^(LSL-Cas9))mice, 2) p53^(flox/flox); T; H11^(LSL-Cas9) (PT; H11^(LSL-Cas9)) mice inwhich virally initiated tumors would lack p53, and 3) Lkb1^(flox/flox);T; H11^(LSL-Cas9) (LT; H11^(LSL-Cas9)) mice in which virally initiatedtumors would lack Lkb1 (FIG. 24a and FIG. 29a ).

LT; H11^(LSL-Cas9) mice were the first to show signs of tumordevelopment including tachypnea and weight loss approximately fivemonths after AAV administration. This is consistent with the rapidgrowth of lung tumors in mice with a Cre-regulated Kras^(G12D) alleleand loss of Lkb1. LT; H11^(LSL-Cas9) mice had very high tumor burdens,resulting from many primary lung tumors (FIG. 24b,c and FIG. 29b-d ).Histological analysis of the lungs of these mice confirmed the presenceof large adenomas and adenocarcinomas (FIG. 24b and FIG. 29 b). PT;H11^(LSL-Cas9) mice also developed numerous large primary lung tumors.Compared to the LT; H11^(LSL-Cas9) mice, tumors initiated in PT;H11^(LSL-Cas9) mice had more pronounced nuclear atypia, a featurecharacteristic of p53-deficiency. Finally, T; H11^(LSL-Cas9) micedeveloped smaller, less advanced lesions, even at later time points(FIG. 24b,c and FIG. 29b-d ). Mice transduced with a 10-fold lower doseof AAV-Kras^(HDR)sg Kras/Cre developed proportionally fewer tumors (FIG.29e ).

Several LT; H11^(LSL-Cas9) and PT; H11^(LSL-Cas9) mice transduced withAAV-Kras^(HDR)/sg Kras/Cre also developed invasive primary lung tumors,disseminated tumor cells (DTCs) in their pleural cavities, and lymphnode metastases (FIG. 24d,e and FIG. 29f,g ). Thus, AAV-Kras^(HDR)/sgKras/Cre-induced tumors can progress into malignant and metastatic lungcancer.

We estimated the efficiency of AAV/Cas9-mediated somatic HDR in the lungby infecting Kras^(LSL-G12D); PT and Kras^(LSL-G12D); LT mice with a1:10,000 dilution of AAV-Kras^(HDR)/sg Kras/Cre, such that oncogenicKras^(G12D) would be expressed in all virally transduced cells. Thesemice developed approximately half as many tumors as mice in whichoncogenic Kras alleles were generated by AAV/Cas9-mediated somatic HDR.This result is consistent with an HDR frequency between 0.02% and 0.1%,enabling the robust initiation of multiple lung tumors in parallel inindividual mice (FIG. 24c ). Importantly, delivery of an analogousvector library without sg Kras (AAV-Kras^(HDR)/Cre) to T, PT, and LTmice did not lead to efficient tumor initiation, suggesting that neitherp53-nor Lkb1-deficiency, combined with high-level AAV vectortransduction, is sufficient to drive lung tumorigenesis (FIG. 24c andFIG. 30).

To verify that tumors initiated using AAV-Kras^(HDR)/sg Kras/Creharbored mutant Kras^(HDR) alleles, we analyzed the Kras locus inFACS-isolated Tomato^(positive) cancer cells from large, individual lungtumors from LT; H11^(LSL-Cas9) and PT; H11^(LSL-Cas9) mice. PCRamplification using primers specific to the Kras^(HDR) allele confirmedthe presence of an oncogenic Kras allele with a unique barcode in eachtumor (FIG. 24f and FIG. 31a,b ). Interestingly, despite the absence ofany detectable HDR bias and the relatively uniform representation ofmutant alleles in the initial AAV library, only five of the thirteenKras variants were identified in ˜50 large lung tumors (FIG. 24f ). Thisresult is consistent with differential selection of Kras variants inlung tumorigenesis.

By analyzing individual tumors, we were able to carefully assess boththe Kras^(HDR) allele as well as the second Kras allele present in tumorcells (FIG. 31). Approximately half of the oncogenic Kras^(HDR) allelesresulted from perfect HDR events, in which a Kras point mutation and aunique barcode were seamlessly recombined into the endogenous Kraslocus. The remaining Kras^(HDR) alleles were seamless from the 5′ endthrough mutant exon 2, but contained small duplications, insertions, ordeletions in intron 2 (FIG. 31d ). Importantly, none of thesealterations would be expected to disrupt splicing from mutant exon 2into exon 3. Additionally, almost all tumors harbored Cas9-inducedindels in the second Kras allele, which is consistent with frequent lossof the wild type KRAS allele in oncogenic KRAS-driven human tumors (FIG.31e,f ). While previous studies have documented enhanced Kras^(G12D)-and Kras^(Q61L)-driven lung tumor growth following inactivation of thewild type Kras allele in mice, our results suggest that many oncogenicKras variants are likely suppressed by wild type Kras during lung tumorgrowth.

In addition to driving human lung cancer, oncogenic KRAS is nearlyubiquitous in human pancreatic ductal adenocarcinoma (PDAC). Expressionof Kras^(G12D) or Kras^(G12V) and inactivation of p53 leads to thedevelopment of PDAC in mouse models. To determine whetherAAV/Cas9-mediated somatic HDR could also induce cancer-initiatingoncogenic point mutations in pancreatic epithelial cells, we transducedPT; H11^(LSL-Cas9) mice with AAV-Kras^(HDR)/sg Kras/Cre by retrogradepancreatic ductal injection (FIG. 25a and FIG. 32a ). These micedeveloped precancerous pancreatic intra-epithelial neoplasias (PanINs)as well as PDAC (FIG. 25b and FIG. 32b,c, f ). Several mice alsodeveloped invasive and metastatic PDAC, consistent with the aggressivenature of the human disease (FIG. 25c and FIG. 32d-f ). Sequencing ofKras^(HDR) alleles from several large pancreatic tumor masses uncoveredoncogenic Kras alleles with unique barcodes (FIG. 24d ). Interestingly,although just four samples were analyzed, only Kras^(G12D) andKras^(G12V) were observed the two most frequent KRAS mutations in humanpancreatic cancer. Consistent with the requirement for oncogenic Kras toinitiate PDAC, transduction of pancreatic cells in PT mice by retrogradepancreatic ductal injection of our negative control AAV-Kras^(HDR)/Crevector did not induce any pancreatic tumors (FIG. 32f ).

Human soft tissue sarcomas also frequently harbor mutations in the RASpathway as well as in TP53. Sarcomas have been induced in geneticallyengineered mice through the expression of Kras^(G12D) and inactivatingp53. To determine whether AAV/Cas9-mediated somatic HDR could be used tointroduce point mutations into Kras and drive sarcoma formation, weperformed intramuscular injections of AAV-Kras^(HDR)/sg Kras/Cre intothe gastrocnemii of PT; H11^(LSL-Cas9) mice (FIG. 25e and FIG. 33a ).These mice developed rapidly growing and invasive sarcomas that harboreduniquely barcoded Kras^(G12D), Kras^(G12A), and KrasG^(13R) alleles(FIG. 25f-h and FIG. 33). The successful application of this platformfor modeling tumorigenesis—from initiation through malignantprogression—in divergent tissues, highlights its broad applicability formultiplexed functional analyses of oncogenic driver mutations in a widerange of cancer types.

Whereas current methods to assess gene function in autochthonous cancermodels largely rely on manual quantification of tumor number and size,we established a simple yet high-throughput and multiplexed approach tolink tumor cell number with tumor genotype directly from bulk tissue(FIGS. 23d and 4a ). Since a unique DNA barcode introduced by HDR into asomatic cell will increase in number as the cell divides, the relativenumber of cancer cells in a given tumor can be determined by deepsequencing of the barcode region. Furthermore, the absolute number ofcells in each tumor can be estimated by adding a normalization controlto each sample prior to deep sequencing. To determine the genotype andestimate the absolute number of cancer cells in each tumor in wholelungs from T; H11^(LSL-Cas9), PT; H11^(LSL-Cas9,) and LT; H11^(LSL-Cas9)mice transduced with AAV-Kras^(HDR)/sg Kras/Cre, we first added DNA from5×10⁵ cells with a known barcode to each sample (FIG. 26a and FIG. 34).We then extracted DNA from the bulk lung samples, PCR-amplified theKras^(HDR) alleles, and deep-sequenced the variant-barcode region ofeach allele (FIGS. 23d and 4a , and FIG. 34).

Following high-throughput sequencing, we corrected for recurrentsequencing errors and the possibility of individual tumors havingidentical barcodes. We then estimated the absolute number of cancercells in each tumor by normalizing tumor barcode sequencing read countsto the number of reads from the normalization control DNA. This analysispipeline was exceptionally reproducible with a high degree ofconcordance in tumor sizes across technical replicates (FIG. 35). Byenabling quantitative analyses of individual tumors in parallel frombulk tissue, this HDR-based barcoding and deep sequencing approachprovides an unprecedented picture of the tumor landscape in vivo.

High-throughput sequencing of the Kras^(HDR) variant-barcode regionuncovered many AAV-Kras^(HDR)/sg Kras/Cre-induced lung tumors in T;H11^(LSL-Cas9), PT; H11^(LSL-Cas9) and LT; H11^(LSL-Cas9) mice (FIG.36a-c ). Normalizing tumor number to the initial representation of eachKras^(HDR) allele in the AAV-Kras^(HDR)/sg Kras/Cre vector libraryallowed us to directly compare the in vivo oncogenicity of each Krasvariant (FIG. 26b and FIG. 36d,e ). Across more than 500 tumors,Kras^(G12D) was the most common variant, consistent with KRAS^(G12D)being the most frequent KRAS mutation in human lung adenocarcinoma innon-smokers. Kras^(G12A), Kras^(G12C), and Kras^(G12V) (the mostfrequent KRAS variants in human lung adenocarcinoma after KRAS^(G12D))as well as KrasG^(13S) were identified as moderate drivers of lungtumorigenesis, but were present in significantly fewer tumors thanKras^(G12D) (FIG. 26b ). Interestingly, Kras^(G12R) and KrasG^(13R) werealso identified as potent oncogenic variants, despite being infrequentlymutated in human lung cancer (FIG. 26b ).

We initiated tumors in PT; H11^(LSL-Cas9), and LT; H11^(LSL-Cas9) miceto directly assess whether concurrent tumor suppressor alterationsmodulate the ability of different Kras variants to initiate and drivetumor growth. Interestingly, although the overall spectrum of Krasoncogenicity changed significantly with Lkb1 inactivation, we did notobserve dramatic differences in the relative tumorigenic potential ofindividual Kras variants in tumors with coincident inactivation of p53or Lkb1 (FIG. 26c-e and FIG. 36). This data is consistent with a modelin which the strength of signaling induced by these oncogenic Krasvariants in vivo is insufficient to engage the p53-pathway; thus, whilep53 functions to constrain tumor progression, it does not limit theinitial expansion of tumors with certain Kras genotypes. Additionally,while Lkb1-deficiency increases tumor growth, the signaling induced byLkb1-deficiency does not preferentially synergize with the downstreamsignals induced by specific mutant forms of Kras.

Since our tumor barcoding and sequencing platform allowed us to identifymany individual lung tumors in parallel from bulk lungs, we anticipatedthat we could also use this approach to overcome the challenge ofidentifying and analyzing individual pancreas tumor clones in multifocaltumor masses initiated in autochthonous mouse models of human PDAC³¹.Therefore, we also analyzed bulk pancreatic tumor samples from PT;H11^(LSL-Cas9) mice transduced with AAV-Kras^(HDR)/sg Kras/Cre (FIG. 26fand FIG. 37a,b ). Barcode sequencing of pancreatic tumor massesuncovered multiple primary tumor clones per mouse, each harboring aKras^(HDR) allele with a point mutation in Kras codon 12 or 13 and aunique DNA barcode. Pancreatic tumors demonstrated oncogenic Kras allelepreferences with Kras^(G12D), Kras^(G12V), and Kras^(G12R) being themost prevalent variants (FIG. 26f ). Notably, these three Kras variantsare also the most prevalent oncogenic KRAS mutations in human PDAC.

In addition to determining the in vivo oncogenicity of specific Krasvariants, our barcode-sequencing approach allowed us to identifycontiguous tumor clones from multi-region sequencing of PDAC masses, anduncover clonal relationships between primary tumors and their metastaticdescendants (FIG. 26g and FIG. 37c ).

The prevalence of a mutation in human cancer is a function of both thefrequency with which the mutation is incurred and the degree to whichthe mutation drives tumorigenesis. By using AAV/Cas9-mediated somaticHDR to introduce point mutations into the endogenous Kras locus in anunbiased manner, we determined that Kras variants have quantitativelydifferent abilities to drive lung tumorigenesis (FIG. 4b and FIG. 36).Furthermore, pancreatic tumors initiated in mice using our HDR-basedapproach demonstrated selection for the same dominant Kras variants ashuman PDACs, suggesting that the spectrum of KRAS mutations observed inhuman PDAC is likely driven by biochemical differences between KRASmutants rather than by differences in their mutation rates (FIG. 26f andFIG. 37).

To begin to understand how the biochemical properties of each Krasvariant influences its in vivo oncogenicity, we investigated therelationship between previously documented biochemical behaviors of Krasvariants and their ability to drive lung or pancreatic tumor formationin our studies (FIG. 38). Notably, although KRAS mutations lead todramatic differences in biochemical features thought to be critical toKRAS function (for example, GTPase activity and RAF kinase affinity), nosingle biochemical property predicted in vivo Kras variant oncogenicity.This result suggests that the in vivo oncogenicity of Kras variants maybe best described by an alternate biochemical property, or perhaps morelikely, through the integration of multiple biochemical outputs.

This work highlights our AAV/Cas9-mediated somatic HDR approach as aquantitative, scalable, and modular approach for cost-effective andsystematic studies of the in vivo oncogenicities of panels of mutationsin parallel. Multiplexed approaches that enable the genetic dissectionof Kras function in vivo represent a critical complement to ongoingbiochemical and cell culture studies of mutant forms of RAS proteins.This method will enable an unprecedented understanding of the functionof diverse mutations in prevalent oncogenes as well as rare, putativelyoncogenic mutations across many common cancer types. Finally, weenvision that this platform will dramatically accelerate both thediscovery and pre-clinical validation of targeted therapies forprecisely defined genetic subtypes of cancer.

Methods

Design, Generation, and Screening of sgRNAs Targeting Kras

To obtain an sgRNA targeting Kras to enhance homology-directed repair(HDR) in somatic mouse cells, we identified all possible 20-bp sgRNAs(using the consensus Cas9 PAM: NGG) targeting Kras exon 2 and theflanking intronic sequences and scored them for predicted on-targetcutting efficiency using an available sgRNA design/scoring algorithm. Wethen empirically determined the cutting efficiency of three sgRNAstargeting Kras (sg Kras#1: GCAGCGTTACCTCTATCGTA; sg Kras#2:GCTAATTCAGAATCACTTTG; sg Kras#3: GACTGAGTATAAACTTGTGG) (FIG. 27a ).Briefly, Lenti-U6-sgRNA/Cre vectors were generated for each sgRNAtargeting Kras as previously described. Q5® site-directed mutagenesis(NEB) was used to insert the sgRNAs into a parental lentiviral vectorcontaining a U6 promoter to drive sgRNA transcription as well as a PGKpromoter driving Cre-recombinase. The cutting efficiency of each sg Kraswas determined via transduction of LSL-YFP; Cas9 cells in culture witheach Lenti-sg Kras/Cre virus. We isolated YFP^(positive) cells 48 hourspost-infection (transduction) by FACS, extracted DNA, PCR-amplified thetargeted Kras locus (forward primer: TCCCCTCTTGGTGCCTGTGTG; reverseprimer: AAGCCCTTCCTGCTAATCTCGGAG) and Sanger-sequenced the amplicons(sequencing primer: GCACGGATGGCATCTTGGACC). Sequencing traces wereanalyzed by TIDE to determine percent indel induction. Since all threesgRNAs induced indels at the anticipated loci, the sg Kras targeting thesequence closest to Kras codons 12 and 13 (sg Kras#3) was used for allsubsequent experiments as this was expected to best facilitate HDR atthe desired locus (FIG. 27a ).

Design, Construction, and Validation of AAV-Kras^(HDR) Plasmid Libraries

Generation of the AAV-Kras^(HDR)/sg Kras/Cre Backbone

The U6-sg Kras/PGK-Cre cassette from pLL3.3; U6-sg Kras/PGK-Cre wasPCR-amplified with Q5® polymerase (NEB), TOPO-cloned (Invitrogen), andverified by sequencing. To generate the AAV-sg Kras/Cre vector, thesequence between the ITRs of the 388-MCS AAV plasmid backbone wasremoved using XhoI/SpeI. The U6-sg Kras/PGK-Cre cassette was digestedfrom the TOPO vector with XhoI/XbaI and the 1.9-kb fragment was ligatedinto the XhoI/SpeI-digested 388-MCS backbone, destroying the SpeI site.A BGH polyA sequence was inserted 3′ of Cre following MluI digestion. A˜2-kb region surrounding exon two of murine Kras was PCR-amplified fromgenomic DNA (forward primer:GCCGCCATGGCAGTTCTTTTGTATCCATTTGTCTCTTTATCTGC; reverse primer:GCCGCTCGAGCTCTTGTGTGTATGAAGACAGTGACACTG). Amplicons were subsequentlycloned into a TOPO vector (Invitrogen). AvrII/BsiWI sites wereintroduced into the TOPO-cloned 2-kb Kras sequence using Q5®site-directed mutagenesis (NEB) (AvrII forward primer:TGAGTGTTAAAATATTGATAAAGTTTTTG; AvrII reverse primer:CCTagGTGTGTAAAACTCTAAGATATTCC; BsiWI forward primer:CTTGTAAAGGACGGCAGCC; BsiWI reverse primer: CGtACGCAGACTGTAGAGCAGC;restriction sites are underlined with mismatching bases in lowercase).The Kras fragment harboring AvrII/BsiWI sites was released from TOPOwith NcoI/XhoI and ligated into NcoI/XhoI-digested AAV-sg Kras/Cre toproduce the AAV-Kras^(HDR)/sg Kras/Cre backbone.

Generation of the AAV-Kras^(HDR)/Cre Backbone

To generate an AAV-Kras^(HDR) backbone without the sgRNA targeting Kras,PGK-Cre was excised from a TOPO clone with NotI/XbaI, and ligated intoNotI/XbaI-digested 388-MCS AAV plasmid backbone. A BGH polyA sequenceand the mouse Kras fragment were added as described above to produce thecontrol AAV-Kras^(HDR)/Cre backbone.

Design and Synthesis of the Diverse Kras Variant/Barcode Region

To introduce a library of activating single point mutations and a DNAbarcode into the Kras^(HDR) sequence of the AAV backbones, wesynthesized four 295-bp Kras fragments with a degenerate “N” base (A, T,C, or G) at each of the first two basepairs of Kras codons 12 and 13(Integrated DNA Technologies) (FIG. 27b ). By design, each of the fourfragment pools consisted of three non-synonymous, single nucleotidemutations at codons 12 and 13 as well as the wild type Kras sequence toserve as a control. Thus, since each of the four pools contained wildtype fragments, the overall representation of wild type Kras alleles wasexpected to be approximately four times higher than each of the mutantKras alleles. The synthesized fragments also contained silent mutationswithin the sg Kras target sequence and the associated protospaceradjacent motif (PAM*), and an eight-nucleotide random barcode created byintroducing degenerate bases into the wobble positions of the downstreamKras codons for individual tumor barcoding (FIG. 27b ). Finally, eachfragment included flanking AvrII and BsiWI restriction sites for cloninginto the AAV-Kras^(HDR) backbones (FIG. 27b ).

Ligation of Kras Mutant/Barcode Fragments into the AAV-Kras^(HDR)Vectors

The four synthesized fragment pools were combined at equal ratios andPCR-amplified (forward primer: CACACCTAGGTGAGTGTTAAAATATTG; reverseprimer: GTAGCTCACTAGTGGTCGCC). Amplicons were digested with AvrII/BsiWI,purified by ethanol precipitation, and ligated into both AAV-Kras^(HDR)backbones (FIG. 27c ). Each ligated plasmid library was transformed intoStbl3 electro-competent cells (NEB) and plated onto 20 LB-Amp plates,which generated ˜3×10⁵ bacterial colonies per library. Colonies werescraped into LB-Amp liquid media and expanded for six hours at 37° C. toincrease plasmid yields to obtain enough plasmid DNA for AAV production.Plasmid DNA was then extracted from bacterial cultures using a Maxiprepkit (Qiagen).

Validation of AAV-Kras^(HDR) Plasmid Libraries

To determine the representation of each Kras variant and thedistribution of barcode nucleotides within each AAV plasm id library,purified AAV plasmid libraries were PCR-amplified with primers tailedwith Illumina adapters (lowercase) containing multiplexing tags(underlined N's) (forward primer:aatgatacggcgaccaccgagatctacactctttccctacacgacgctcttccgatctCTGCTGAAAATGACTGAGTATAAACTAGTAGTC; reverse primer:caagcagaagacggcatacgagatNNNNNNgtgactggagttcagacgtgtgctcttccgatcCTGCCGTCCITTACAAGCGTACG), and then deep-sequenced on a MiSeq (Illumina®).

AAV Capsid Serotype Lung Epithelial Cells Transduction Analysis

Recombinant AAV-GFP vectors were produced using a Ca₃(PO₄)₂ tripletransfection protocol with pAd5 helper, ssAAV-RSV-GFP transfer vectorand pseudotyping plasm ids for each of nine capsids of interest: AAV1,2, 3b, 4, 5, 6, 8, 9_hu14 and DJ. Viruses were produced in HEK293T cells(ATCC) followed by double cesium chloride density gradient purificationand dialysis as previously described. rAAV vector preparations weretitered by TaqMan qPCR for GFP (forward primer: GACGTAAACGGCCACAAGTT;reverse primers: GAACTTCAGGGTCAGCTTGC; probe:6-FAM/CGAGGGCGATGCCACCTACG/BHQ-1). To identify an optimal AAV serotypefor adult lung epithelial cell transduction, each mouse received 60 μlof pseudotyped AAV-GFP at maximal titer via intratrachealadministration. Mice were analyzed 5 days after AAV administration.Lungs were dissociated into single-cell suspensions and prepared forFACS analysis of GFP^(positive) cells as described previously.GFP^(positive) percent ages were determined by analyzing >10,000live-gated cells (see FIG. 28).

Production and titering of AAV-Kras^(HDR) plasmid libraries

AAV libraries were produced using a Ca₃(PO₄)₂ triple transfectionprotocol with pAd5 helper, pAAV2/8 packaging plasmid and the barcodedKras library transfer vector pools described above. Transfections wereperformed in HEK293T cells followed by double cesium chloride densitygradient purification and dialysis as previously described. AAVlibraries were titered by TaqMan qPCR for Cre (forward primer:TTTGTTGCCCTTTATTGCAG; reverse primer: CCCTTGCGGTATTCTTTGTT; probe:6-FAM/TGCAGTTGTTGGCTCCAACAC/BHQ-1).

In Vitro AAV/Cas9-Mediated HDR

The nucleotide changes surrounding the mutations at codon 12 and 13(three nucleotide changes 5′ of codons 12/13 to mutate the sgRNArecognition site and PAM motif, and up to 10 changes in the barcodesequence) made it unlikely that the point mutations at Kras codons 12and 13 would differentially affect the rate of HDR. We neverthelesstested whether HDR efficiency might be influenced by differences in thesequence of individual Kras^(HDR) alleles. To induce in vitroAAV/Cas9-mediated HDR, we transduced LSL-YFP/Cas9 cells with thepurified AAV-Kras^(HDR)/sg Kras/Cre library (FIG. 27e ). Cells weremaintained in cell culture media with 10 μM SCR7 (Xcessbio), aninhibitor of non-homologous end joining (NHEJ), to promote homologydirected repair. 96 hours after transduction, DNA was isolated from theLSL-YFP/Cas9 cells by phenol/chloroform extraction followed by ethanolprecipitation. The Kras locus was amplified from this DNA using a PCRstrategy that we developed for the specific amplification of Kras^(HDR)alleles integrated into the endogenous Kras locus. We thendeep-sequenced these amplicons to determine the representation ofKras^(HDR) alleles following in vitro HDR (see the “Illumina® librarypreparation and sequencing of tumor barcodes from bulk tissue” sectionbelow for details on PCR and sequencing).

Mice and Tumor Initiation

Lkb1^(flox) (L), p53^(flox) (P), R26^(LSL-Tomato) (T) H11^(LSL-Cas9),and Kras^(LSL-G12D) (K) mice have been previously described. AAVadministration by intratracheal inhalation to initiate lung tumors,retrograde pancreatic ductal injection to initiate pancreatic tumors,and intramuscular gastrocnemius injection to initiate sarcomas wasperformed as described. Lung tumors were initiated in PT;H11^(LSL-Cas9), LT; H11^(LSL-Cas9) and T; H11^(LSL-Cas9) mice with 60 μlof AAV-Kras^(HDR)/sg Kras/Cre (1.4×10¹² vg/ml), in PT, LT, and T micewith 60 μl of AAV-Kras^(HDR)/Cre (2.4×10¹² vg/ml), or in KPT and KLTmice with 60 μl AAV-Kras^(HDR)/sg Kras/Cre (1.4×10¹² vg/ml) diluted1:10,000 in 1×PBS. Pancreatic tumors were initiated in PT;H11^(LSL-Cas9) mice with 100-150 μl of AAV-Kras^(HDR)/sg Kras/Cre(1.4×10¹² vg/ml) or in PT mice with 100-150 μl of AAV-Kras^(HDR)/Cre(2.4×10¹² vg/ml). A 1:10 dilution of AAV-Kras^(HDR)/sg Kras/Cre in 1×PBSwas also administered to the lungs or pancreata of mice where indicated.Sarcomas were initiated in PT; H11^(LSL-Cas9) with 30 μl ofAAV-Kras^(HDR)/sg Kras/Cre (5.2×10¹² vg/ml). Mice were euthanized whenthey displayed symptoms of tumor development. The Institutional AnimalCare & Use Committee of Stanford University approved all mouseprocedures.

Analysis of Individual Tumors

Analysis of Individual Lung Tumors

Lung tumor-bearing mice displaying symptoms of tumor development andwere analyzed 4-10 months after viral administration. Lung tumor burdenwas assessed by lung weight and by quantification of macroscopicTomato^(positive) tumors under a fluorescence dissecting scope asindicated (a single LT; H11^(LSL-Cas9) mouse had minimalTomato^(positive) signal that was restricted to a small region of onelung lobe, indicative of improper intratracheal administration of AAV,and was removed from the study). The largest individual lung tumors thatwere not visibly multifocal were dissected from bulk lungs under afluorescence dissecting microscope for sequencing. For some lung tumors,the Tomato^(positive) tumor cells were purified using FACS machines(Aria sorter; BD Biosciences) within the Stanford Shared FACS Facility.Several lung lobes from individual mice were also collected forhistological analysis.

Analysis of Individual Pancreatic Tumor Masses

Pancreatic tumor-bearing mice displayed symptoms of tumor developmentand were analyzed 3-4 months after viral administration. Sincepancreatic tumors largely appeared to be multifocal, individual regionsof the pancreas containing Tomato^(positive) tumor masses were dissectedand FACS-purified for sequencing (a mouse treated with a 1:10 dilutionof AAV-Kras^(HDR)/sg Kras/Cre library also developed pancreatic tumormasses and therefore was included in these analyses). Regions of severalpancreata were kept for histological analysis.

Analysis of Individual Sarcomas

Sarcoma-bearing mice with obvious tumor development and were analyzed3-7 months after viral administration. A region of each sarcoma was keptfor sequencing and an adjacent region was saved for histologicalanalysis.

Characterization of Kras Alleles in Individual Tumors

DNA for sequencing was extracted from FACS-purified tumor cells andunsorted tumor samples with a DNeasy Blood and Tissue Extraction kit(Qiagen). To identify Kras point mutations and barcodes in tumors, wePCR-amplified and sequenced the Kras^(HDR) alleles using two protocolsthat were optimized across several variables including annealingtemperature, extension time, and primer sequences (Protocol 1 forwardprimer: CTGCTGAAAATGACTGAGTATAAACTAGTAGTC; reverse primer:AGCAGTTGGCCTTTAATTGGTT; sequencing primer:AATGATACGGCGACCACCGAGATCTACAC; annealing temperature 66° C.; Protocol2—forward primer: GCTGAAAATGACTGAGTATAAACTAGTAGTC; reverse primer:TTAGCAGTTGGCCTTTAATTGG; sequencing primer: GCACGGATGGCATCTTGGACC;annealing temperature: 64° C.). These protocols were used tospecifically amplify integrated Kras^(HDR) alleles from individualtumors as each incorporated a forward primer overlapping the engineeredmutations in the PAM region upstream of codons 12 and 13, and a reverseprimer outside the homology arm. Long extension times (2-3 minutes) wereused to enable amplification of all Kras^(HDR) alleles, even thosecontaining insertions or duplications in intron 2 of the Kras locus(FIG. 31d ).

Apart from introducing the desired point mutations into the endogenousKras locus via HDR, targeting Kras exon 2 using CRISPR/Cas9 was alsoexpected to result in indels at the cut site following DNA repair byNHEJ instead of HDR. To characterize these modifications, we used ageneric PCR protocol to amplify both Kras alleles (forward primer:TCCCCTCTTGGTGCCTGTGTG; reverse primer: GGCTGGCTGCCGTCCTTTAC; sequencingprimer: CAAGCTCATGCGGGTGTGTC; annealing temperature: 72° C.). A spectrumof insertions and deletions at the site of DNA cleavage was identifiedby this approach (FIG. 31e,f ). For some individual tumor samples, thesequence of both Kras alleles was not immediately obvious following theabove PCR and sequencing strategies. PCR products from these sampleswere TOPO cloned (Invitrogen) and transformed, and several colonies fromeach sample were plasmid prepped and sequenced to characterize both Krasalleles in each tumor. This approach was reproducible and reliableacross both biological and technical replicates, and a variety ofHDR-induced oncogenic Kras alleles were identified. Indel-containingKras alleles were identified in ˜50 tumors (FIG. 31a,b ).

These analyses also uncovered several other unexpected features in someof the Kras alleles from individual lung tumors. Three distinct missensemutations at codon 24 (I24L, I24N, I24M) were observed in a small subsetof the individual lung tumors analyzed. The function of thesealterations, if any, is unknown.

Furthermore, we initially anticipated that recombination of theKras^(HDR) template into the endogenous Kras locus would occur outsideof the AvrII and BsiWI sites engineered into the Kras^(HDR) template(FIG. 31c ). However, the AvrII site, engineered by altering 2 basepairs 97 base pairs upstream of exon 2, was absent in 5 out of 25 tumorsin which we directly analyzed this region of the Kras^(HDR) allele. TheBsiWI site, engineered by altering 1 base pair 20 base pairs downstreamof exon 2, was absent in 11 out of 58 tumors. These findings indicatedthat while recombination of the Kras^(HDR) template most often occurredwithin the larger, more distal homology arms, it also occurred at adetectable frequency within very short regions of homology that areflanked by 5′ and 3′ mismatches (including the PAM* mutation, a Krascodons 12 or 13 mutation, and 8 potential mismatches within thebarcode).

After we initially identified the presence of duplications in theKras^(HDR) alleles in some tumors, we designed PCR primers tospecifically amplify duplications of Kras exon 2 that occurred on eitherside of the HDR-integrated Kras locus (Right-hand duplication—forwardprimer: TGACCCTACGATAGAGGTAACG; reverse primer: CTCATCCACAAAGTGATTCTGA;sequencing primer: TGACCCTACGATAGAGGTAACG; Left-hand duplication—forwardprimer: TGAGTGTTAAAATATTGATAAAGTTTTTG; reverse primer:TCCGAATTCAGTGACTACAGATG; sequencing primer:TGAGTGTTAAAATATTGATAAAGTTTTTG). Each of these duplication-specific PCRprotocols used adjacent primer pairs in opposite orientations, ensuringthat amplification would only occur if a duplication was present.Duplications of varying lengths were identified (FIG. 31d ), includingduplications of the second half of wild type exon 2 or the entire exon 2(but lacking critical regions of the splice acceptor). Deletions andduplications of regions of intron 2 were also observed. Furthermore, weobserved integrations of parts of the AAV vector, including the U6promoter and viral ITR, into intron 2. Given the size and location ofthese alterations, none would be expected to change splicing of Krasmutant exon 2 to exon 3, consistent with the requirement of oncogenicKras to drive tumorigenesis.

Generating a Normalization Control for High-Throughput Sequencing from aCell Line with a Known Kras^(HDR) Allele and Barcode

To establish a cell line to use as a sequencing normalization control, asingle large tumor was dissected from an PT; H11^(LSL-Cas9) mouse,digested into a single cell suspension, and plated to generate a cellline. After expanding these cells and then extracting DNA, Kras exon 2was PCR amplified (forward primer: TCCCCTCTTGGTGCCTGTGTG; reverseprimer: GGCTGGCTGCCGTCCTTTAC). The PCR product was sequenced (usingspecific and generic sequencing primers described above) to confirm thepresence of a Kras^(HDR) allele and a barcode. A single Kras^(G12V)allele with a unique barcode (CGGGAAGTCGGCGCTTACGATC) was identified.The genomic DNA from this cell lines was used as a normalization controlfor high-throughput sequencing for all bulk lung samples (FIG. 34).

Bulk Tissue Processing and DNA Extraction

Bulk Lung Tissue Processing and DNA Extraction

Bulk lung samples were dissected from infected (transdcued) mice andstored at −80° C. prior to processing. To extract DNA for sequencing,samples were thawed and transferred to 50 mL conical tubes. 20 mLs oflysis buffer (100 mM NaCl, 20 mM Tris pH7.6, 10 mM EDTA pH8.0, 0.5% SDSin H₂O) and 200 μL Proteinase K (20 mg/mL) were added to each sample.Next, 3 μg (˜5×10⁵ genomes) of normalization control DNA was added toeach sample (FIG. 25 a and FIG. 35b ). Samples were then carefullyhomogenized using a tissue blender, which was cleaned between eachsample by progressing through clean 10% bleach, 70% ethanol, and 1×PBS.Homogenized samples were lysed at 55° C. overnight. DNA was isolatedfrom tissue lysates by phenol/chloroform extraction followed by ethanolprecipitation (FIG. 37a,b ).

Bulk Pancreatic Tissue Processing and DNA Extraction

Pancreatic tumor masses were dissected, digested, and viable(DAPI^(negative)), lineage (CD45, CD31, Ter119, F4/80)^(negative),Tomato^(positive) cells were isolated by FACS. No normalization controlwas added to the pancreatic cancer samples. DNA was isolated from theFACS isolated neoplastic cells using a DNeasy Blood and TissueExtraction kit (Qiagen), and then further purified by ethanolprecipitation.

Illumina® Library Preparation and Sequencing of Tumor Barcodes from BulkTissue

To uncover the number and size of tumors harboring each Kras variant ina massively parallel and quantitative manner, we developed a two-roundPCR strategy that enabled multiplexed Illumina® sequencing of barcodedKras^(HDR) alleles (FIG. 27f ). For the 1^(st) round of PCR, we used aforward primer complementary to the Kras^(HDR) sequence containing thethree PAM and sgRNA target site mutations (PAM*; bold in the 1^(st)round forward primer sequence) (1^(st) round forward primer:GCTGAAAATGACTGAGTATAAACTAGTAGTC) (SEQ ID NO: 2), and a reverse a primercomplementary to a downstream region of the endogenous Kras locus notpresent in the HDR template in the AAV-Kras^(HDR)/sg Kras/Cre vector(1^(st) round reverse primer: TTAGCAGTTGGCCTTTAATTGG) (SEQ ID NO: 3).This primer pair was chosen to specifically amplify genomic Kras^(HDR)alleles without amplifying abundant wild type Kras alleles or potentialepisomal AAV-Kras^(HDR)/sg Kras/Cre vectors present in DNA purified frombulk tumor-bearing tissue. Additionally, a P5 adapter (italicized), 8-bpcustom i5 index (N's), and Illumina® sequencing primer sequence (read 1)(underlined) was included at the 5′ end of the 1^(st) round forwardprimer to enable multiplexed Illumina® sequencing (1^(st) round forwardprimer for Illumina sequencing:

(SEQ ID NO: 4) AATGATACGGCGACCACCGAGATCTACACNNNNNNNNACACTCTTTCCCTACACGACGCTCTTCCGATCTGCTGAAAA TGACTGAGTATAAACTAGTAGTC).

Importantly, since the characterization of Kras^(HDR) alleles inindividual tumors uncovered some variability in HDR resulting in diverseindels in Kras intron 2 (FIG. 32d ), only 4 (lung samples) or 6(pancreas samples) cycles were performed in the 1^(st) round of PCR tominimize the potential for bias during the amplification of products ofvariable length. Furthermore, a high-efficiency polymerase (Q5® HotStart High-Fidelity polymerase, NEB; 64° C. annealing temperature) and along extension time (3:00 minutes) were used to ensure robustamplification of all Kras^(HDR) alleles. Kras^(HDR) alleles in genomiclung DNA were amplified using between 4 and 40 separate 100 μL PCRreactions and then pooled following amplification to reduce the effectsof PCR jackpotting (FIG. 34a ). Each of these 100-μL PCR reactionscontained 4 μg of DNA template to amplify from a large initial pool ofKras^(HDR) alleles. Following the 1^(st) round of amplification, allreplicate PCR reactions were pooled and 100 μL of each sample wascleaned up using a QIAquick PCR Purification Kit (Qiagen).

Purified 1^(st) round PCR amplicons were used as template DNA for a 100μL 2^(nd) round Illumina® library PCR (Q5® Hot Start High-Fidelitypolymerase, NEB; 72° C. annealing temperature; 35 cycles for lungsamples, 40 cycles for pancreas samples). The 2^(nd) round of PCRamplified a 112-bp region entirely within the Kras exon 2 sequencepresent in 1^(st) round PCR amplicons. The 2^(nd) round reverse primercontained a P7 adapter (italicized), reverse complemented 8-bp custom i7index (“Ns”), and reverse complemented Illumina sequencing primersequence (read 2) (underlined) at the 5′ end to enable dual-indexed,paired-end sequencing of Illumina libraries (2^(nd) round reverse primer#1: CAAGCAGAAGACGGCATACGAGATNNNNNNNNGTGACTGGACTTCAGACGTGTGCTCTTCCGATCCGTAGGGTCATACTCATCCACA) (SEQ ID NO: 5). The 2^(nd) round PCRforward primer was complementary to the P5 Illumina adapter added to theamplified Kras^(HDR) allele by the forward primer during the 1^(st)round PCR (2^(nd) round forward primer: AATGATACGGCGACCACCGAGATCTACAC)(SEQ ID NO: 6). This primer was used to amplify 1^(st) round PCRamplicons without amplifying any contaminating genomic DNA that may havebeen carried over from the 1^(st) round PCR reaction. Furthermore, asecond reverse primer encoding the P7 adaptor sequence was added to the2^(nd) round PCR reaction at the same concentration as the two otherprimers (2^(nd) round reverse primer #2: CAAGCAGAAGACGGCATACGAGAT) (SEQID NO: 7). This primer binds the reverse complemented P7 adaptorsequence added to the Kras^(HDR) amplicons by 2^(nd) round reverseprimer #1. Since the 2^(nd) round PCR was performed over 35-40 cycles,the P7 adaptor (2^(nd) round reverse primer #2) was added to limit theamount of non-specific amplification produced by the lengthy 2^(nd)round reverse primer #1.

After the 2^(nd) round of amplification, 100-μL PCR reactions were runon a 2.5% agarose gel and a band of the expected size was excised. DNAwas extracted from gel fragments using a QIAquick Gel Extraction Kit(Qiagen). The quality and concentration of the purified Illumina®libraries was determined using a Bioanalyzer (Agilent). IndividualIllumina® libraries with unique dual-indices were then pooled togethersuch that libraries originally derived from mice with greater tumorburden were represented at a higher ratio in the final pool than thosefrom mice with lower tumor burdens (FIG. 34a ). A total of 35 individualsamples were combined into two Illumina® library pools. The quality andconcentration of each pool was confirmed on a Bioanalyzer (Agilent).Each final Illumina® library pool was then deep-sequenced on anIllumina® HiSeq lane using multiplexed, 150 bp paired-end Rapid Runsequencing program (Elim Biopharmaceuticals).

Analysis of Illumina Sequencing Data to Estimate the Size and Number ofBarcoded Tumors

We developed a pipeline to call tumors from our de-multiplexed Illumina®sequencing data. The pipeline tallies unique barcode sequences andeliminates recurrent sequencing errors using an algorithm designed todenoise deep-sequencing data of amplicons (DADA2). We tailored thisalgorithm to minimize the occurrence of spurious tumor calls, andminimize technical biases (including variation in read depth, variationin Illumina® sequencing machine error rates, and variation in barcodediversity). This pipeline, including modifications for the analysis oftumor genotypes and barcodes following AAV/Cas9-mediated somaticHDR-driven tumorigenesis, is described below.

Merging, Filtering, and Trimming Paired-End Reads

Although our Illumina® sequencing libraries contained a small 112-bpfragment of the Kras^(HDR) alleles, we performed 150 bp paired-endsequencing of these fragments and merged the overlapping forward andreverse reads to reduce the likelihood of Illumina® sequencing errors inKras codons 12 and 13 and the barcode region of the Kras^(HDR) alleles.Overlapping paired end-reads were merged, quality-filtered, and trimmedusing PANDAseq (fragment length: 60 bp; forward trimming primer:ATGACTGAGTATAAACT; reverse trimming primer: CTCATCCACAAAGTGA).

Calling Unique Tumors

Even after merging forward and reverse reads to reduce sequencingerrors, an average of ˜1 error per 10,000 bases was detected, presumablyfrom recurrent Illumina® sequencing errors (or less likely fromrecurrent PCR errors). Given this error rate, we expected that readsfrom a large, uniquely barcoded tumor containing single nucleotidemismatches would be called as small, spurious tumors of ˜ 1/10,000^(th)the size of the large real tumor. This phenomenon was discernible byeye, as we observed small clusters of spurious “tumors” that were ˜3-4orders of magnitude smaller and contained 1 nucleotide deviationsrelative to the largest tumor in specific mice. Additionally, eachKras^(HDR) variant-barcode pair also generated recurrent sequencingerrors in the mutant base in the oncogenic codon 12 or 13.

To accurately call tumors, we developed a computational and statisticalpipeline for the analysis of tumor barcode sequencing data with thefollowing steps:

Training an Error Model from Non-Barcode Regions of Reads and ClusteringUnique Read Pileups into Tumors Using DADA2

We estimated the residual rate of sequencing/PCR error from the 7nucleotides upstream of KRAS codon 12 and the 7 nucleotides downstreamof the final barcode base. We then used our model of sequencing errorsto cluster unique read pileups (truncated to within 7 nucleotides of thebarcoded bases) into unique tumors via DADA2. A minimum confidence inunique origin of the clusters of 0.01 (i.e. omega_a=0.01) was used. Alarger threshold increased the number of unique tumors called in a mousesample. We chose this larger value because paired-end sequencingappeared to give us greater confidence that unique read pileups weretruly distinct tumors. For example, we found that this thresholdeliminated all unintended read sequences (e.g. reads with inappropriatenucleotides outside of the barcode), and that this threshold called atotal number of lesions within each mouse that was more consistentbetween biological replicates. These were important considerations sincewithout proper handling of read errors the number of called tumors canpositively correlate with sequencing read depth. Finally, we removed anytumors with DNA sequences that deviated by only 1 nucleotide from alesion that was 10,000× larger. This affected only 1.56% of tumor calls.

Normalizing Read Pileups to a Normalization Control to Get ApproximateTumor Size

After generating the read pileups and performing the correctionsdescribed above, we normalized the number of reads from each calledtumor to the number of reads from the normalization control that wasspiked into each sample prior to DNA extraction from bulk tumor-bearingtissue lysates. This allowed us to generate a reasonable estimate of thenumber of cells in each tumor and allowed us to merge data from mice ofthe same genotype and treatment. However, there are several factors thatimpact our ability to accurately quantify the absolute number of cellswithin each tumor.

A first consideration is that some of the Kras^(HDR) alleles inindividual tumors harbored insertions or deletion in Kras intron 2,inside the PCR primers for Illumina® sequencing. Although the presenceof different sized amplicons could generate a PCR bias, we attempted toreduce this by performing only 4-6 cycles in the 1^(st) round ofIllumina library PCR, using a long extension time (˜3 minutes), andusing a fast (20-30 seconds/kb), high fidelity polymerase (OS®; NEB). Asthe final Illumina® library PCR product in 2^(nd) round of amplificationis short and uniform across all samples, PCR amplication should not bebiased in this step.

Furthermore, given that the Kras variants and barcodes are knocked intothe endogenous Kras locus, it is possible that in some tumors thisregion is genomically amplified (which has been documented inKras^(G12D)-driven lung tumors initiated in mouse models of lungcancer). Although Kras amplifications do not typically result in veryhigh Kras copy number, any amplification would lead to a slightoverestimation of the number of cells in tumors with amplifiedKras^(HDR) alleles since our conversion from read count to cell numberassumes that each cell contains a single copy of the barcoded Kras^(HDR)allele.

Lastly, the normalization control itself was generated from cells from atumor with a known duplication in Kras intron 2, which produces a largerPCR product in the 1^(st) round of the Illumina library preparation thantumors without a duplication. Thus, any PCR bias away from the Krasalleles in the normalization control would result in a systematicunderestimation of the size of tumors without duplications.

Estimating the Barcode Overlap Rate and Correcting Tumor SizeDistributions

Sequencing of tumor barcodes from 35 samples on two lanes of anIllumina® Hi-Seq Rapid Run, combined with our analysis pipeline, enabledthe detection of unique barcodes with read counts covering over fiveorders of magnitude. Thus, the unprecedented resolution of this approachenables the detection of large lesions as well as small hyperplasiaswithin bulk tissue. However, the ability to detect a large numbers oflesions within bulk tissue increases the probability of barcodecollisions: the occurrence of two or more lesions with the same DNAbarcode in the same mouse. Barcode collisions can overstate the size ofobserved tumors because two small “colliding” tumors would be identifiedas a single, larger tumor. Therefore, we developed a statistical modelof barcode collisions to ensure that this issue was modest and did notovertly bias the estimated sizes of called tumors.

Our model of barcode collisions accounts for the likelihood p ofobserving each of the 24,576 possible barcodes i for each Kras variantin our study. A majority of the reproducible variation between barcodefrequencies in our pool derives from statistically independent variationin the nucleotide frequencies at each wobble base (i.e. each barcode isnot equally likely in the pool because there was subtle variation innucleotide concentrations during synthesis of the barcodes fragments)(FIG. 23g ). Thus, we estimate the independent frequency f_(b,n) of eachnucleotide n at every base b in the barcode and use this table topredict barcode likelihoods based on each barcode's sequence B_(i, b, n)(where B is 1 if barcode i possesses nucleotide n at position b and 0otherwise) as follows:

$p_{i} = {\underset{b}{\Pi}{B_{i,b,n} \cdot f_{b}^{n}}}$

Here, matrix notation is used to denote a dot product. This modelpredicts every barcode's frequency with only 21 free parameters. Becausesome residual over-representation of barcodes persisted in the lungsamples, we simply discarded the 10% most frequently observed barcodes,after correcting for nucleotide frequencies, from all lung analyses.These most frequently observed barcodes were identified independent ofour mouse experiments by Illumina® sequencing (MiSeq) of ourAAV-Kras^(HDR)/sg Kras/Cre plasmid pool prior to virus production. Afterthis processing, we then renormalized Σ_(i)p_(i) to one.

We then assumed that the occurrences of each barcode within each mousewas a Multinomial sampling process. The mean number of collisions C_(i)for each observed barcode within each mouse is then:

$\begin{matrix}{{C_{i}\left( {p_{i},N} \right)} = {\sum\limits_{k = 1}^{\infty}\; {\left( {k - 1} \right){P\left( {{k;p_{i}},N} \right)}}}} \\{= {\mu_{i} - {P\left( {{{k = 0};p_{i}},N} \right)} - 1}} \\{= {{Np}_{i} + \left( {1 - p_{i}} \right)^{N} - 1}}\end{matrix}$

Here, μ_(i) denotes the mean number of barcodes within each mouse, whileN denotes the total number of tumors (both unknowns). N is determinedfrom the observed number of tumors in each mouse N^((obs)) using theequation N^((obs))=N−Σ_(i)C_(i)(p_(i), N) and Brent's Method.

This model found that barcode collisions were generally rare in ourmouse samples (on average 4.04%). However, the likelihood of collisionscan vary by mouse and by Kras variant. For example, the averagepredicted number of collisions for WT Kras^(HDR) alleles was 5.8% and ashigh as 12% in one mouse. WT Kras^(HDR) alleles were expected toexperience the highest number of collisions since WT Kras vectors wereintentionally represented ˜4 fold more than each mutant Kras vector inthe initial AAV-Kras^(HDR)/sg Kras/Cre plasmid pool) (FIG. 23f ). Thus,we divided the size of each lesion by 1+C_(i) to minimize the bias thatbarcode collisions impart on tumor size distributions. Becausecollisions are rare events, the particular number of collisions within aparticular mouse can differ substantially from C_(i). Because of thislimitation, we believe that this correction minimizes systematic bias intumor size distributions resulting from barcode collisions; however, itcannot effectively identify the specific collisions that occurred.

Determining Illumina® Sequencing Quality and Reproducibility

To determine whether Kras variants had quantitatively differentabilities to drive tumorigenesis, we elected to focus on tumorsestimated to contain more than 100,000 cells (i.e. ⅕^(th) the “size” ofthe normalization control DNA added to each sample that was derived from˜5×10⁵ cells). Regression analysis of tumors above this cell numbercutoff from replicate samples (independent sample preparation,sequencing, and processing) demonstrated high correlation (all R² valueswere above 0.99; see FIG. 36). Furthermore, the estimated number ofcells in tumors below this cutoff were more likely to be biased bybarcode collisions and variability in PCR amplification and sequencing,all of which would have decreased our ability to accurately call thesize and number of tumors harboring each Kras variant.

Analysis of Sequencing Data from Bulk Tumor-Bearing Lungs

We quantified the relative number of tumors harboring each Kras variantby counting tumors above 100,000 cells in all mouse genotypes with aH11^(LSL-Cas9) allele (PT; H11^(LSL-Cas9), LT; H11^(LSL-Cas9), and T;H11^(LSL-Cas9)), and dividing each variant by its initial representationin the AAV-Kras^(HDR)/sg Kras/Cre plasmid pool (for this analysis, theinitial representation of each variant in the plasmid pool wascalculated from the total number of reads associated with each Krasvariant after removing barcodes above the 98^(th) percentile of barcodeabundance; this restriction did not appreciably alter results, and wassimply applied to ensure that extremely abundant variant-barcode pairsdid not overtly impact the overall representation of specific variants).

Relative tumor number was then scaled such that WT Kras variants had arepresentation of 1. There was a relatively small number of WTKras^(HDR) alleles that appeared to arise from tumors above 100,000cells. These could represent tumors in which an HDR event created thenon-oncogenic Kras WT genotype but which nonetheless evolved into atumor for other reasons, or the WT Kras variant ‘hitchhikes’ with anoncogenic Kras variant by co-incident HDR in the same lung cell followedbut expansion driven by the oncogenic variant.

A small number of residual cells from individual tumors that weredissected from bulk tissue (and analyzed as described above) wereusually detectable in our bulk tumor sequencing data. In all analyses oftumor size these dissected tumors were excluded, as we could not infertheir true size. However, when analyzing the number of tumors above100,000 cells in each treated mouse genotype, we included data fromindividually dissected tumors since dissectible tumors were always amongthe largest observed within any mouse and, therefore, certainly abovethe 100,000 cell threshold.

Statistically significant differences in tumor number were determinedusing Fischer's Exact Test. For each variant two tests were performed,comparing to either the frequency of G12D or WT Kras^(HDR) alleles. Allp-values are Bonferroni-Corrected for the number of variantsinvestigated and are two-sided. A two-sided “many cells” Pearson's ChiSquared Test was used to compare the distribution of tumor numbersacross all Kras variants in PT; H₁₁ ^(LSL-Cas9) and LT; H₁₁ ^(LSL-Cas9)mice relative to T; H11^(LSL-Cas9) mice.

Example 3. The Fitness Landscape of Tumor Suppression in LungAdenocarcinoma in Vivo

The functional impact of most genomic alterations found in cancer, aloneor in combination, remains largely unknown. With experiments describedherein, integration of tumor barcoding, CRISPR/Cas9-mediated genomeediting, and ultra-deep barcode sequencing is demonstrated forinterrogating pairwise combinations of tumor suppressor alterations inautochthonous mouse models of human lung adenocarcinoma. The tumorsuppressive effects of 31 common lung adenocarcinoma genotypes aremapped, revealing a rugged landscape of context-dependence anddifferential effect strengths.

Results

Cancer growth is largely the consequence of multiple, cooperativegenomic alterations. Cancer genome sequencing has catalogued a multitudeof alterations within human cancers, however the combinatorial effectsof these alterations on tumor growth is largely unknown. Most putativedrivers are altered in less than ten percent of tumors, suggesting thatthese alterations may be inert, weakly-beneficial, or beneficial only incertain genomic contexts. Inferring genetic interactions throughco-occurrence rates alone is practically impossible, as the number ofpossible combinations scales factorially with candidate gene number.Genetically engineered mouse models can provide insight into genefunction in tumors growing within an autochthonous setting, howeverpractical considerations have prevented broad studies of combinatorialtumor suppressor gene inactivation (FIG. 41). Hence, our understandingof the genetic interactions that drive tumor growth in vivo remainslimited.

To address these practical challenges, a method was developed (describedherein) to quantitatively measure the effect of many different tumorsuppressor gene alterations in parallel using tumor barcoding coupledwith deep-sequencing (Tuba-seq). Tuba-seq combines geneticallyengineered mouse models of lung adenocarcinoma with tumor suppressorinactivation (e.g., CRISPR/Cas9-mediated), tumor barcoding, anddeep-sequencing. Because Tuba-seq measures the size of every tumor andis compatible with multiplexing tumor genotypes in individual mice,growth effects can be measured with unprecedented precision,sensitivity, and throughput. Here, this approach is employed The growthof oncogenic Kras^(G12D)-driven lung tumors with 31 common tumorsuppressor genotypes is quantified (FIG. 39). Unexpected geneticinteractions were identified, the effects of most tumor suppressors werefound to be context-dependent, and several patterns of geneticalterations in human lung adenocarcinomas were explained.

The tumor suppressor TP53 is inactivated in more than half of human lungadenocarcinomas. To determine the effect of p53 deletion on the growthsuppressive effects of ten other putative tumor suppressors, tumors wereinitiated in Kras^(LSL-G12D); Rosa26^(LSL-tdTomato); H11^(LSL-Cas9) (KT;Cas9) and KT; p53^(flox/flox); Cas9 (KPT; Cas9) mice using a pool ofbarcoded Lenti-sg RNA/Cre vectors targeting many common tumor suppressorgenes and four barcoded Lenti-sg Inert/Cre vectors (Lenti-sgTS-Pool/Cre; FIGS. 39. 41, and 42). Barcodes contained two componentsthat uniquely identify each tumor and its sgRNA (sgID-BC; FIG. 42). Thenumber of neoplastic cells in each tumor of each genotype was determined15 weeks after tumor initiation when the lungs contained widespreadhyperplasias, adenomas, and some early adenocarcinomas. The sgID-BCregion was amplified from bulk tumor-bearing lung genomic DNA, theproduct was deep sequenced, and the Tuba-seq analysis pipeline(described herein) was applied.

Tuba-seq analysis of KT; Cas9 and KPT; Cas9 mice uncovered an alteredspectrum of tumor suppressive effects for many of the genes in oursurvey (FIGS. 39 and 43). Tumor sizes were summarized by twopreviously-vetted measures: Lognormal (LN) mean and the size of the 95thpercentile tumor (FIG. 39). In p53-deficient tumors, inactivation ofRb1, Setd2, Lkb1/Stk11, Cdkn2a, or Apc still provided a growthadvantage, while Smad4, Arid1a, and Atm emerged as tumor suppressorsonly in the absence of p53 (FIGS. 39 and 43). The emergence ofadditional tumor suppressors in this background suggests thatp53-deficiency potentiates subsequent tumor evolution. By allowing moremutations to be adaptive, p53 loss may decrease the predictability oftumor evolution and facilitate future tumor evolution, including theemergence of treatment resistance and metastatic disease.

Coincident deletion of p53 not only allowed more alterations to beadaptive, but also significantly changed the magnitude of effect oftumor suppressor loss. In KT; Cas9 mice, Rb1-deficiency increased tumorsize less than Lkb1- or Setd2-deficiency (FIG. 39 and FIG. 44a ;P<0.0001 bootstrap test unless otherwise specified). In contrast, in thep53-deficient background, Rb1-deficiency conferred a growth advantagecomparable to that of Lkb1- or Setd2-deficiency (P>0.05), consistentwith a strong complementary interaction between the p53 and Rb1 tumorsuppressor pathways (FIG. 39). Quantification of Cas9-generated indelsat each targeted locus in bulk KPT; Cas9 lung DNA confirmed comparablyhigh percentages of Lkb1, Setd2, and Rb1 alleles with indels (FIGS. 39and 45). Finally, the effect of co-incident inactivation of p53 and Rb1on lung cancer growth was confirmed using conventional Cre/LoxP-basedmouse models (FIG. 44).

The quantitatively-different growth benefits of Rb1 inactivation in p53proficient versus p53 deficient tumors presented the opportunity toinvestigate whether changes in the fitness strength of a driver altersthe frequency of its alterations in human lung adenocarcinomas. Indeed,co-occurrence of RB1 alterations (SNVs and CNVs) and TP53 alterationswere enriched in human lung adenocarcinoma (P=0.03; FIG. 39 and FIG.44). Notably, despite a ˜5-fold enrichment in the co-occurrence of thesetwo alterations, this interaction would be statistically insignificantin a nave survey of all potential pairwise driver interactions aftercorrecting for multiple-hypothesis testing, thus illustrating the needto functionally study genetic interactions beyond co-occurrence patterns(P=0.32 after Bonferroni correction for 10 pairwise interactions).

Next, the effects of combinatorial loss of Lkb1 and other putative tumorsuppressors was investigated by initiating tumors withLenti-sgTS-Pool/Cre in KT; Lkb1^(flox/flox); Cas9 (KLT; Cas9) mice(FIGS. 40 and 43). Lkb1 was investigated because it dramaticallyincreases lung tumor growth in autochthonous models and is frequentlyinactivated in human lung adenocarcinoma (FIG. 41). Interestingly, boththe number of adaptive tumor suppressor losses and the median growthbenefit was attenuated in the already fast-growing Lkb1-deficient tumors(irrespective of changes in statistical power between mouse backgrounds,P<0.05, Methods). This once again demonstrates that a single alterationcan change the entire fitness landscape of tumors. General attenuationof fitness benefits, termed diminishing returns epistasis, is common inevolution, and suggests that tumors may eventually reach a fitnessplateau.

Apc and Rb1 inactivation were the only alterations that provided asignificant growth advantage to Lkb1-deficient tumors (FIG. 40). Theability of Rb1-deficiency to increase tumor size, even with coincidentLkb1-deficiency, emphasizes the integral role of Rb1 in cell cycleregulation and fundamentally different mechanism of action from Lkb1loss. Apc loss is also a key driver of lung cancer growth, and Apc wastumor suppressive in all three backgrounds studied.

Surprisingly, the effect of Setd2-deficiency on the growth ofLkb1-deficient tumors was modest and statistically insignificant (FIG.40). This redundancy was striking because both Lkb1 and Setd2inactivation strongly promote growth in KT; Cas9 and KPT; Cas9 mice andbecause there is no evidence that these genes function in the samepathway. Thus, the context dependence of Setd2 inactivation was testedand confirmed by initiating tumors with Lenti-sgNeo2/Cre and Lenti-sgSetd2/Cre in KPT, KPT; Cas9, and KLT; Cas9 mice. Setd2 inactivationenhanced Lkb1-proficient lung tumor growth, while conferring little, ifany, growth advantage to Lkb1-deficient tumors (P<0.05 of (sg Setd2 inKPT; Cas9/KLT; Cas9)/(sgNeo2 in KPT; Cas9/KLT; Cas9) for histologicalanalysis, P<0.0001 for Tuba-seq analysis, FIGS. 40, and 46). Thisobservation is also well supported by the mutual exclusivity ofLKB1/STK11 and SETD2 alterations in human lung adenocarcinoma (P<0.001,FIGS. 40 and 46).

Most genes in these studies exhibited context-dependent growth effects,driving tumor growth only in the presence or absence of p53 or Lkb1(FIG. 40). Even the tumor suppressor alterations that conferredadvantage in all three contexts (Rb1 and Apc) still exhibitedcontext-dependent magnitudes of tumor suppression. Such wide-spreadcontext-dependency is overlooked by global surveys of drivers, wheredriver interactions are either ignored or presumed to be sufficientlyrare and/or weak to justify considering only marginal correlations.Nonetheless, our fitness measurements overall agree with mutationco-occurrence patterns in human lung cancer, despite the limitedstatistical resolution of these data (Spearman R=0.50, P<0.05, FIG. 47).Furthermore, while lung cancers do not appear to be unique in theirdegree of context-dependency (FIG. 47) and the findings here suggestthat direct measurement of context-dependency in other cancer types iswarranted.

This rugged landscape of tumor evolution has several implications.First, to understand gene function, it can be important to investigateputative drivers in multiple genetic contexts, as most genes in thesurvey (8 of 11) were only adaptive in some contexts (FIG. 40). Second,broader fitness profiling is desirable. The power analyses here suggestthat ˜500 moderate-strength interactions could be surveyed usingTuba-seq with a one hundred-mouse cohort (FIG. 48). Larger genomicscreens could survey more putative drivers, interactions with otheroncogenic events, multiple sgRNAs targeting the same gene, or tripletsof tumor suppressor alterations. Lastly, this extensivecontext-dependency suggests that most driver alterations sweep tofixation infrequently because they are beneficial only in specificgenetic contexts.

The studies described herein of the fitness effects of combinatorialtumor suppressor losses in vivo identified unexpected geneticinteractions that were validated by traditional methodologies as well asby human lung adenocarcinoma genomics data. The barcoded and multiplexedgenome-editing approach described herein could easily be utilized tointerrogate the functional consequences of these genetic interactions,including their impact on therapeutic response, cell signaling, and/ormetastatic progression.

1-31. (canceled)
 32. A method of measuring tumor size for a plurality ofclonally independent tumors of the same tissue, the method comprising:(a) contacting a tissue with a plurality of barcoded nucleic acid cellmarkers, thereby generating a plurality of distinguishable lineages ofheritably marked neoplastic cells within the contacted tissue; (b) aftersufficient time has passed for at least a portion of the heritablymarked neoplastic cells to undergo at least one round of division,performing high-throughput nucleic acid sequencing to detect and measurequantities of at least two of the of barcoded nucleic acid cell markerspresent in the contacted tissue, thereby generating a set of measuredvalues; and (c) calculating, using the set of measured values as input,a number of heritably marked neoplastic cells present in the contactedtissue for at least two of said distinguishable lineages of heritablymarked neoplastic cells.
 33. The method of claim 32, wherein said tissuecomprises neoplastic cells and/or tumors prior to step (a).
 34. Themethod of claim 32, wherein the high-throughput nucleic acid sequencingof step (b) is performed on a biological sample collected from thetissue.
 35. The method of claim 32, wherein the high-throughput nucleicacid sequencing of step (b) is performed on a tissue sample of thecontacted tissue.
 36. The method of claim 32, wherein each barcodednucleic acid cell marker of the plurality of barcoded nucleic acid cellmarkers corresponds to a known cell genotype for a lineage of heritablymarked neoplastic cells.
 37. The method of claim 32, wherein saidcontacting comprises genetically altering cells of the tissue togenerate the heritably marked neoplastic cells.
 38. The method of claim32, wherein the barcodcd nucleic acids induce neoplastic cell formation.39. The method of claim 32, wherein the barcodcd nucleic acids induceneoplastic cell formation and include one or more of: homology directedrepair (HDR) DNA donor templates, nucleic acids encoding one or moreoncogenes, nucleic acids encoding one or more wildtype proteins, nucleicacids encoding one or more mutant proteins, nucleic acids encodingCRISPR/Cas guide RNAs, nucleic acids encoding short hairpin RNAs(shRNAs), or nucleic acids encoding a genome editing protein.
 40. Themethod of claim 39, wherein the genome editing protein is selected from:a CRISPR/Cas RNA-guided protein, a CRISPR/Cas RNA-guided protein fusedto a transcriptional activator or repressor polypeptide, a Cas9 protein,a Cas9 protein fused to a transcriptional activator or repressorpolypeptide, a zinc finger nuclease (ZFN), a TALEN, a phage-derivedintegrase, a Cre protein, a Flp protein, and a meganuclease protein. 41.The method of claim 32, wherein the barcodcd nucleic acids are linear orcircular DNA molecules.
 42. The method of claim 32, wherein the barcodcdnucleic acids are selected from: plasmids, synthesized nucleic acidfragments, and minicircles.
 43. The method of claim 32, wherein thebarcodcd nucleic acids are RNA/DNA hybrids or nucleic acid/proteincomplexes.
 44. The method of claim 32, wherein the tissue is aninvertebrate tissue.
 45. The method of claim 32, wherein the tissue is avertebrate tissue.
 46. The method of claim 32, wherein the tissue is amammalian or a fish tissue.
 47. The method of claim 32, wherein thetissue is a rat tissue, a mouse tissue, a pig tissue, a non-humanprimate tissue, or a human tissue.
 48. The method of claim 32, whereinthe tissue is part of a living animal.
 49. The method of claim 32,wherein the tissue is an engineered tissue grown outside of an animal.50. The method of claim 32, wherein the tissue is selected from: muscle,lung, bronchus, pancreas, breast, liver, bile duct, gallbladder, kidney,spleen, blood, gut, brain, bone, bladder, prostate, ovary, eye, nose,tongue, mouth, pharynx, larynx, thyroid, fat, esophagus, stomach, smallintestine, colon, rectum, adrenal gland, soft tissue, smooth muscle,vasculature, cartilage, lymphatics, prostate, heart, skin, retina,reproductive system, and genital system.
 51. The method of claim 32,wherein after sufficient time has passed for at least a portion of theheritably marked neoplastic cells to undergo at least one round ofdivision, the method further comprises: (i) detecting and/or measuring abiomarker of the heritably marked neoplastic cells, and (ii)categorizing the heritably marked neoplastic cells based on the resultsof said detecting and/or measuring of the biomarker.
 52. The method ofclaim 51, wherein the biomarker one or more of: cell proliferationstatus, cell type, developmental cell lineage, cell death, or cellularsignaling state.
 53. The method of claim 32, wherein the cell marker isdelivered to the tissue via viral vector.
 54. The method of claim 53,wherein the viral vector is selected from: a lentiviral vector, anadenoviral vector, an adeno-associated viral (AAV) vector, a bocavirusvector, a foamy virus vector, and a retroviral vector.